Dynamics of Fragment Formation in Neutron Rich Matter

TL;DR

This study investigates the formation and dynamics of nuclear fragments in neutron-rich matter during expansion, simulating conditions similar to neutron star mergers, revealing new phases and fragmentation behaviors.

Contribution

It introduces a novel phase called pregnocchi and analyzes fragment formation dynamics in expanding neutron-rich matter using molecular dynamics simulations.

Findings

Reproduced known pasta phases and identified a new pregnocchi phase.

Discovered different fragmentation regimes based on initial conditions.

Observed neutron clouds around clusters influencing r-process configurations.

Abstract

Background: Neutron stars are astronomical systems with nucleons submitted to extreme conditions. Due to the longer range coulomb repulsion between protons, the system has structural inhomogeneities. Several interactions tailored to reproduce nuclear matter plus screened Coulomb term reproduce these inhomogeneities known as nuclear pasta. These structural inhomogeneities, located in the crust of neutron stars, can also arise in expanding systems depending on the thermodynamic conditions (temperature, proton fraction, ...) and the expansion velocity. Purpose: We aim to find the dynamics of the fragments formation for expanding systems simulated according to the little big bang model. This expansion resembles the evolution of neutron stars merger. Method: We study the dynamics of the nucleons with semiclassical molecular dynamics models. Starting with an equilibrium configuration, we expand the system homogeneously until we arrive to an asymptotic configuration (i. e. very low final densities). We study, with four different cluster recognition algorithms, the fragment distribution throughout this expansion and the dynamics of the cluster formation. Results: Studying the topology of the equilibrium states, before the expansion, we reproduced the known pasta phases plus a novel phase we called pregnocchi, consisting of proton aggregates embedded in a neutron sea. We have identified different fragmentation regimes, depending on the initial temperature and fragment velocity. In particular, for the already mentioned pregnocchi, a neutron cloud surrounds the clusters during the early stages of the expansion, resulting in systems that give rise to configurations compatibles with the emergence of r-proccess. Conclusions: These calculations pave the way to a comparision between Earth experiments and neutron star studies.

Tables

Table 1. Table 1 : Three different expansion of Neutron Star Matter: Lasagna (fast expansion): x = 0.4 𝑥 0.4 x=0.4 , η ˙ = 0.01 ​ s − 1 ˙ 𝜂 0.01 superscript s 1 \dot{\eta}=0.01\text{s}^{-1} , T = 0.8 ​ MeV 𝑇 0.8 MeV T=0.8\,\text{MeV} ; Lasagna (slow expansion): x = 0.4 𝑥 0.4 x=0.4 , η ˙ = 0.0001 ​ s − 1 ˙ 𝜂 0.0001 superscript s 1 \dot{\eta}=0.0001\text{s}^{-1} , T = 0.8 ​ MeV 𝑇 0.8 MeV T=0.8\,\text{MeV} ; Pregnocchi: x = 0.1 𝑥 0.1 x=0.1 , η ˙ = 0.0001 ​ s − 1 ˙ 𝜂 0.0001 superscript s 1 \dot{\eta}=0.0001\text{s}^{-1} , T = 0.1 ​ MeV 𝑇 0.1 MeV T=0.1\,\text{MeV}

	Lasagna (fast expansion)	Lasagna (slow expansion)	Pregnocchi
$\rho =0.05{\text{fm}}^{-3}$
$\rho =0.0001{\text{fm}}^{-3}$
$\rho =0.00003{\text{fm}}^{-3}$

Equations

$$V_{np}^{\text{CMD}}(r)$$

$$V_{nn}^{\text{CMD}}(r)$$

$$V_{TF}(r)=q^2\frac{e^{-r/\lambda }}{r}$$

$$i\in C_n^{\text{MST}}\iff \exists j\in C_n\mid r_{ij}<r_{cut}$$

$$i\in C_n^{\text{MSTE}}\iff \exists j\in C_n^{\text{MSTE}}:V_{ij}+K_{ij}\le 0$$

$${ϵ}_n=\sum _{i\in C_n}{K}_i^{CM}+\sum _{i,j\in C_n}{V}_{ij}$$

$$L(t)=L_0(1+\dot{\eta }t)$$

$$\mathbf{v}={\mathbf{v}}_0+\dot{\eta }{\mathbf{r}}_0$$

$$V_{ij}(r)=\{\begin{array}{ll}-V_0& \text{if }r\le a\\ 0& \text{if }r>a.\end{array}$$

$$K_{\text{CM}}$$

$$=\frac{n^3}{3}m{v}^2+\mathcal{O}(n^2)$$

$$V_{\text{CM}}$$

$$=-2n^2{V}_0$$

Full text

Neutron stars are astronomical systems with nucleons submitted to extreme conditions. Due to the longer range coulomb repulsion between protons, the system has structural inhomogeneities. Several interactions tailored to reproduce nuclear matter plus screened Coulomb term reproduce these inhomogeneities known as nuclear pasta. These structural inhomogeneities, located in the crust of neutron stars, can also arise in expanding systems depending on the thermodynamic conditions (temperature, proton fraction, …) and the expansion velocity.

Purpose

We aim to find the dynamics of the fragments formation for expanding systems simulated according to the little big bang model. This expansion resembles the evolution of neutron stars merger.

Method

We study the dynamics of the nucleons with semiclassical molecular dynamics models. Starting with an equilibrium configuration, we expand the system homogeneously until we arrive to an asymptotic configuration (i. e. very low final densities). We study, with four different cluster recognition algorithms, the fragment distribution throughout this expansion and the dynamics of the cluster formation.

Results

Studying the topology of the equilibrium states, before the expansion, we reproduced the known pasta phases plus a novel phase we called pregnocchi, consisting of proton aggregates embedded in a neutron sea. We have identified different fragmentation regimes, depending on the initial temperature and fragment velocity. In particular, for the already mentioned pregnocchi, a neutron cloud surrounds the clusters during the early stages of the expansion, resulting in systems that give rise to configurations compatibles with the emergence of r-proccess.

Conclusions

These calculations pave the way to a comparision between Earth experiments and neutron star studies.

I Introduction

Neutron rich matter is present in several astronomical objects in the universe, for example: neutron stars, proto-neutron stars and core-collapse supernovae. The supernova explosion of a massive star, combined with gravitational collapse, compresses the core up to densities of atomic nuclei. This gives rise to a system known as proto-neutron star, which eventually ends up in a neutron star.

The neutron rich environment also gives rise to the possibility of a rapid neutron capture, r-process, that consists of the rapid capture of neutrons. The r-process is fundamental to understand the abundancy of heavy elements, and several places have been candidates for it to happen. Supernovae have been prime candidates for long, but recent observations and models hint that also neutron star mergers can yield r-process Eichler et al. (1989); Freiburghaus et al. (1999); Tanvir et al. (2013). The compression of neutron star matter as a possible source for r-process nuclei was first discussed in Ref. Lattimer and Schramm (1974). According to hydrodynamic models Goriely et al. (2011), these have typically velocity gradients of $\dot{\eta}=10^{-21}\,\text{c/fm}<\dot{\eta}<4\cdot 10^{-20}\,\text{c/fm}$.

The original works of Ravenhall et al. Ravenhall et al. (1983) and Hashimoto et al. Hashimoto et al. (1984) used a compressible liquid drop model to study neutron rich matter, and have shown that the now known as the pasta phases –lasagna, spaghetti and gnocchi– are solutions to the ground state of neutron star matter. The study of neutron rich matter has since been approached with different models, which show that nuclear pasta arises due to the interplay between nuclear and Coulomb forces in an infinite medium. We classify the different approaches in two large groups: mean field and microscopic.

Mean field works include the Liquid Drop Model, by Lattimer et al. Page et al. (2004), Thomas-Fermi, by Williams and Koonin Williams and Koonin (1985), among others Oyamatsu (1993); Lorenz et al. (1993); Cheng et al. (1997); Watanabe et al. (2000); Watanabe and Iida (2003); Nakazato et al. (2009). Microscopic models include Quantum Molecular Dynamics, used by Maruyama et al. Maruyama et al. (1998); Kido et al. (2000) and by Watanabe *et al.*Watanabe et al. (2003), Simple Semiclassical Potential, by Horowitz et al. Horowitz et al. (2004a) and Classical Molecular Dynamics, used in our previous works Dorso et al. (2012).

In some recent studies, phases different from the typical nuclear pasta were found. The work by Nakazato et al. Nakazato et al. (2009), inspired by polymer systems, found also gyroid and double-diamond structures, with a compressible liquid drop model. Dorso et al. Dorso et al. (2012) obtained pasta phases different from those already mentioned with molecular dynamics, studying mostly their characterization at very low temperatures. In our previous work Alcain et al. (2014a) we have shown that these new pasta phases had an opacity peak (i. e., a local maximum in the opacity) in the characteristic wavelength of the Urca neutrinos for symmetrical neutron star matter.

Among the advantages of classical and semiclassical models are the accessibility to position and momentum of all particles at all times, which allows the calculation of correlations of all orders. Moreover, no specific structure is hardcoded in the model, as it happens with most mean field models. This enables the study of the structure of the nuclear medium from a particle-wise point of view. Many models exist with this goal, including quantum molecular dynamics Maruyama et al. (1998), simple-semiclassical potential Horowitz et al. (2004a) and classical molecular dynamics Lenk et al. (1990). In these models the Pauli repulsion between nucleons of equal isospin is hard-coded in the interaction. On the other hand, a specific Pauli potential developed in Dorso et al. (1987) was used in the QCNM Dorso and Randrup (1988) and later in Ref. Hartnack et al. (1989).

The relative inaccessibility to these astronomical objects means a restriction in the observables available. One of them, studied extensively in the recent years, is the neutrino opacity and the mean free path Horowitz et al. (2004b, a); Alcain and Dorso (2017). In this work, we study another possible observable from the neutron rich matter: the result of the fragmentation of neutron rich matter, related to the already mentione r-process. Multifragmentation in nuclear systems has been studied before Bonasera et al. (2000); Chikazumi et al. (2001), but mostly with nuclear matter (without Coulomb interaction). In a recent work by Caplan et al Caplan et al. (2015), expanding neutron star matter has been studied as possible explanations for nucleosynthesis in neutron star mergers.

In Section II we introduce the model used along this work, that includes the potential parametrization (II.1) and the Coulomb interaction (II.2). Section III describes the different cluster recognition algorithms used in this work, and Section IV explains how we simulate the expansion of the system. Finally, we draw conclusions in Section VI. In the appendix we perform a detailed analysis on the stability of one of the cluster recognition algorithms.

II The Model

II.1 Classical Molecular Dynamics

In this work, we study fragmentation of Neutron Star Matter under pasta-like conditions with model similar to the classical molecular dynamics model CMD. CMD has been used in several heavy-ion reaction studies to: help understand experimental data Chernomoretz et al. (2002); identify phase-transition signals and other critical phenomena López and Dorso (2000); Barranon et al. (2001); Dorso and López (2001); Barranón et al. (2003); Barrañón et al. (2007); and explore the caloric curve Barrañón et al. (2004) and isoscaling Dorso et al. (2006, 2011). CMD uses two two-body potentials to describe the interaction of nucleons, which are a combination of Yukawa potentials:

[TABLE]

where $V_{np}$ is the potential between a neutron and a proton, and $V_{nn}$ is the repulsive interaction between either $nn$ or $pp$. The cutoff radius is $r_{c}=5.4\,\text{fm}$ and for $r>r_{c}$ both potentials are set to zero. The Yukawa parameters $\mu_{r}$, $\mu_{a}$ and $\mu_{0}$ were determined to yield an equilibrium density of $\rho_{0}=0.16\,\text{fm}^{-3}$, a binding energy $E(\rho_{0})=16\,\text{MeV/nucleon}$ and a compressibility of $250\,\text{MeV}$.

Based on this model, we developed a new set of parameters that yield the same values for $\rho_{0}$, $E(\rho_{0})$ and compressibility, which we called New Medium. We show in figure 1 an example that compares the proton-neutron potential for the different models and the developed for this work: SSP, CMD and New Medium.

To simulate an infinite medium, we used this potential with $N=5500$ particles under periodic boundary conditions, with different proton fractions (i. e. with $x=Z/A=0.2<x<0.4$) in cubical boxes with sizes adjusted to have densities $\rho=0.05\,\text{fm}^{-3}$ and $\rho=0.08\,\text{fm}^{-3}$. These simulations have been done with LAMMPS Plimpton (1995), using its GPU package Brown et al. (2012).

II.1.1 Ground State Nuclei

Although the $T=0$ state of this classical nuclear matter at normal densities is a simple cubic solid, nuclear systems can be mimicked by adding enough kinetic energy to the nucleons. To study nuclei, for instance, liquid-like spherical drops with the right number of protons and neutrons are constructed confined in a steep spherical potential and then brought to the ground state by cooling them slowly from a rather high temperature until they reach a self-contained state. Removing the confining potential, the system is further cooled down until a reasonable binding energy is attained. The remaining kinetic energy of the nucleons helps to resemble the Fermi motion. To compare the different microscopic models used throughout the literature with the New Medium, we show in figure 2 the binding energies of ground-state nuclei obtained with CMD, SSP and New Medium.

II.2 Coulomb interaction in the model

Since a neutralizing electron gas embeds the nucleons in the neutron star crust, the Coulomb forces among protons are screened. We model this screening effect with the Thomas-Fermi approximation, used with various nuclear models Maruyama et al. (1998); Dorso et al. (2012); Horowitz et al. (2004b). According to this approximation, protons interact via a Yukawa-like potential, with a screening length $\lambda$:

[TABLE]

Theoretical estimates for the screening length $\lambda$ are $\lambda\sim 100\,\text{fm}$ Fetter and Walecka (2003), but we set the screening length to $\lambda=20\,\text{fm}$. This choice was based on previous studies Alcain et al. (2014b), where we have shown that this value is enough to adequately reproduce the expected length scale of density fluctuations for this model, while larger screening lengths would be a computational difficulty. We analyze the opacity to neutrinos of the structures for different proton fractions and densities.

III Cluster recognition

In typical configurations we have not only the structure known as nuclear pasta, but also a nucleon gas that surrounds the nuclear pasta. In order to properly characterize the pasta phases, we must know which particles belong to the pasta phases and which belong to this gas. To do so, we have to find the clusters that are formed along the simulation.

One of the algorithms to identify cluster formation is Minimum Spanning Tree (MST). In MST algorithm, two particles belong to the same cluster $\{C^{\text{MST}}_{n}\}$ if the relative distance of the particles is less than a cutoff distance $r_{cut}$:

[TABLE]

Based on MST algorithm, and taking into account that typically the neutron rich matter structure is set by the proton backbone, we developed MSTpC: an algorithm that calculates the MST cluster of protons alone and finds the cloud of neutrons that lay within $r_{c}$ of each proton cluster. The MST cluster definition works correctly for systems with no kinetic energy, and it is based in the attractive tail of the nuclear interaction. However, if the particles have non-zero relative momenta, we can have a situation of two particles that are closer than the cutoff radius, but with a large relative kinetic energy.

The problem of fragment recognition in nuclear physics has had a strong development in the last years, specially related to the analysis of numerical simulations of intermediate energy heavy ion collisions. In this case, it is quite clear that the asymptotic state of the system is a very dilute set of fragments with a collective expansion mode and composed by cold fragments. In the asymptotic state, fragments will be far away from each other, and as such, application of the above mentioned MST algorithm yields an accurate description of the fragmentation. However, if one is interested in the analysis of the time evolution of the fragment structure, it is clear that the MST will not provide information because during the expansion nearby particles may have very different momenta, which are not considered in MST definition. Therefore, two particles can be very close to each other –within range of the attractive potential– but with high relative momentum being recognized as a bound pair according to MST. This unwanted behavior can be partially solved using the MSTE algorithm, in which two particles belong to the same cluster $\{C^{\text{MSTE}}_{n}\}$ if they are energy bound:

[TABLE]

This is an approximate solution to our problem, but it shows signs of instability even in some simple cases (see appendix).

One of the most sophisticated methods to find the energetically bound clusters is the Early Cluster Recognition Algorithm (ECRA) Dorso and Randrup (1993). In this algorithm, the particles are partitioned in different disjoint clusters $C^{\text{ECRA}}_{n}$, with the total energy in each cluster:

[TABLE]

where $K^{CM}_{i}$ is the kinetic energy relative to the center of mass of the cluster. The set of clusters $\{C^{\text{ECRA}}_{n}\}$ then is the one that minimizes the sum of all the cluster energies $E_{\text{partition}}=\sum_{n}\epsilon_{n}$.

As mentioned above, expanding systems have the property that the asymptotic state is easy to calculate, so the efficiency of other recognition algorithms apart from MST can be easily tested. Moreover, the quality of the recognition algorithm can be weighted considering how early in the evolution of the expanding system it is able to identify the clusters corresponding to the asymptotic state. In this sense, the ECRA algorithm has shown that it is able to recognize fragments very early in the evolution, providing a new view of the dynamics of the fragment formation i. e., fragments are early formed in exploding systems. ECRA algorithm can be easily used for small systems Dorso and Balonga (1994), but being a combinatorial optimization, it cannot be used in large systems.

To find approximate solutions, the originally proposed method is similar to simulated annealing Dorso and Randrup (1993). Another other choice was developed by Puente Puente (1999), and it introduces a Binary Fusion Model. In this model, the initial configuration is with all clusters being monomers (one particle per cluster). With this starting point, $E_{\text{partition}}^{0}=0$ the steps that follow are:

Explore all potential mergers of two clusters and bookkeep the resulting $E_{\text{partition}}^{i+1}$ from each potential merger. 2. 2.

Pick the merger that results in the lowest $E_{\text{partition}}^{i+1}$. 3. 3.

If $E_{\text{partition}}^{i+1}<E_{\text{partition}}^{i}$, perform the merge and go back to step 1; otherwise, stop iteration.

All of these algorithms for cluster recognition should give the same results for the asymptotic state.

III.1 Infinite Clusters

In Ref. Alcain and Dorso (2016) we developed an algorithm for the recognition of infinite clusters across the boundaries. We explain here in detail the implementation for MST clusters in 2D, being the MSTE and 3D extension straightforward. In figure 3 we see a schematical representation of 2D clusters recognized in a periodic cell, labeled from 1 to 6 (note that these clusters don’t connect yet through the periodic walls).

In order to find the connections of these clusters through the boundaries, we draw a labeled graph of the clusters, where we connect clusters depending on whether they connect or not through a wall and label such connection with the wall label. For example, we begin with cluster 1. It connects with cluster 2 going out through the E wall, therefore we add a $1\rightarrow 2$ connection labeled as E. Symmetrically, we add a $2\rightarrow 1$ connection labeled as W. Now we go for the pair 1–3. It connects going out through the S wall, so we add $1\rightarrow 3$ labeled as S and $3\rightarrow 1$ labeled as N. Cluster 1 does not connect with 4, 5, or 6, therefore those are the only connections we have. Once we’ve done that, we get the graph of figure 4.

We now wonder whether these subgraphs represent an infinite cluster or not. In order to have an infinite clusters, we need to have a loop (the opposite is not true: having a loop is not enough to have an infinite cluster, as we can see in subragph 5–6), so we first identify loops and mark them as candidates for infinite clusters. Every connection adds to a loop (since the graph connections are back and forth), but we know from inspecting the figure 4 that the cluster 1–2–3 is infinite. Finding out what makes, in the graph, the cluster 1–2–3 infinite is key to identify infinite clusters. And the key feature of cluster 1–2–3 is that its loop 1–2–3–1 can be transversed through the walls E–E–S, while loops like 5–6 can be transversed only through E–W. Now, in order for the cluster to be infinite, we need it to extend infinitely in (at least) one direction. So once we have the list of walls of the loop, we create a magnitude $I$ associated to each loop that is created as follows: beginning with $I=0$, we add a value $M_{i}$ if there is (at least one) $i$ wall. The values are: $M_{E}=1$, $M_{W}=-1$, $M_{N}=2$, $M_{S}=-2$. If $I$ is nonzero, then the loop is infinite. For example, for the loop E–E–S, we have E and S walls, so $I=M_{E}+M_{S}=3$ and the loop is infinite. For the loop E–W, $I=M_{E}+M_{W}=0$, and the loop is finite.

IV Expansion

In order to expand the neutron rich matter that simulates an infinite system with periodic boundary conditions, we follow the microscopic big bang method, as explained by Dorso and Strachan in Ref. Dorso and Strachan (1996). It consists of an expansion of the simulation box at a constant isotropic rate:

[TABLE]

where $L$ is the length of the simulation box in every direction and $L_{0}$ is the initial length. With only this box resizing, the system would expand dynamically. To simulate an expansion, we need to also give the particles an extra radial velocity that maches that of the box in the edges of the simulation:

[TABLE]

Since we are working with periodic boundary conditions, when a particle crosses a boundary, we must take into account the original expansion, so we do not only change the particle position but also the velocity. For example, if the particle crosses the left-hand boundary of the periodic box, the velocity of the image particle $v_{i}^{\dagger}$ on the right-hand must be modified $v_{i}^{\dagger}=v_{i}+L_{0}\,\dot{\eta}$. This prescription for an expansion is mathematically equivalent to Hubble’s law in astrophysics Chikazumi et al. (2001).

V Results

V.1 Configuration dependence with the potential

Different models for the interaction yield different equations of state and, consequently, different configurations. For comparison, we show in figure 5 different snapshots for the three models we studied: CMD medium, New medium and SSP. These snapshots are near ground-states, with very low temperature ($T=0.1\,\text{MeV}$), density $\rho=0.05\,\text{fm}^{-3}$ and a proton fraction of $x=0.1$. The differences are very noticeable: while CMD medium potential has no identifiable structure, the New medium and SSP potentials clearly show agglomerations of proton (due to the binding interaction with neutrons) embedded in a neutron sea. This structure is what we shall call pre-gnocchi. This is the first time such a structure has been identified and it is also a very interesting qualitative difference observed among parametrizations of the equation of state.

To compare the potentials in a different configuration, we show in figure 6 different snapshots for the three models we studied: CMD medium, New medium and SSP. These snapshots are near ground-states, with very low temperature ($T=0.1\,\text{MeV}$), density $\rho=0.05\,\text{fm}^{-3}$ and a proton fraction of $x=0.4$.

V.2 Asymptotic Mass Distribution

When the system expands, the structure breaks down into finite fragments. For long enough times, these fragments remain stable (since they don’t interact with each other). We will refer to this as the asymptotic fragments.

We expanded several initial configurations to find their asymptotic mass distribution. For the first example, we show in figure 7 the asymptotic mass distribution (calculated with the MSTE algorithm) for $x=0.1$, $\rho=0.05\,\text{fm}^{-3}$ and $T=0.8\,\text{MeV}$ for two expansion velocities: fast ($\dot{\eta}=0.01\,\text{s}^{-1}$) and slow ($\dot{\eta}=0.0001\,\text{s}^{-1}$). We can see here that the slow expansion allows the existence of fragments with mass of up to 60 (20 of which are protons) while the fast expansion produces smaller fragments of up to 20 (6 protons). This is an expected behavior, since the faster expansion, the larger the excitation energy. Therefore, a faster expansion is supposed to break clusters that would otherwise be stable. A similar behavior can be seen in figure 8, where we expand the system for $x=0.4$, $\rho=0.05\,\text{fm}^{-3}$ and $T=0.1\,\text{MeV}$ for the same fast and slow expansion velocities. Another relevant characteristic of the asymptotic mass distribution –not shown in the figures due to scale limitations– is that the fast expansion has a non negligible fraction of lone neutrons (about $4\%$), while the slow expansion hardly presents any ($0.1\%$).

V.3 Fragment Formation

We now turn to the analysis of some examples of the system evolution in time: when and how are these fragments formed. We take first the expansion of the system with $x=0.4$, $\rho=0.05\,\text{fm}^{-3}$ and $T=0.5\,\text{MeV}$. We show in the first two columns of table 1 the initial and the asymptotic state with the slow and the fast expansion. While the initial condition is an infinite cluster, in the asymptotic regime we have a fragment distribution with many finite clusters. It is interesting to note that the fast expansion resembles a mechanical fracture, in which the fragments are formed within each sheet of the lasagna, while the slow expansion looks more like a thermal expansion in which the asymptotic system loses any resemblance to the original structure. The clusters break into many fragments because their large size cannot withstand the energy associated with the expansion of the system.

A very interesting scenario is the expansion of the system with low proton fraction: $x=0.1$, $\rho=0.05\,\text{fm}^{-3}$ and $T=0.1\,\text{MeV}$. In the third column of table 1 we show both the initial condition and the asymptotic configuration for $\dot{\eta}=0.0001\,\text{s}^{-1}$.

Unlike the previous scenario, there is a clear proton backbone the clusters already exists immersed in a neutron sea. It can be visually identified when we draw the protons with a much larger size than neutrons, as is in this set of figures. As the system expands, it is modified. This raises the question: does the cluster distribution change substantially? The answer to this question requires a deep analysis of the time evolution of the cluster distribution, and we no longer can rely on a visual inspection; we need to use the cluster recognition algorithms. Such an analysis has been performed for finite systems for example in Ref. Dorso and Balonga (1994); Strachan and Dorso (1997). In figure 9 we show the initial and final configuration with the MSTE algorithm. Note that the cluster distribution changes radically in both aspects: the size and the proton fraction. The proton fraction change is to be expected, since as the system expands, less neutrons are within the range of the proton cluster. However, this effect alone does not explain the change of size: while the initial condition shows a cluster of up to 80 protons, the asymptotic condition’s largest cluster is of about 30 protons. Did a cluster break down while the system was expanding?

To analyze this, we study the MST distribution of protons alone, shown in figure 10. According to this figure, we see that the cluster distribution of protons did not change substantially (only one proton cluster broke down) and effectively the largest cluster has 32 protons. Does the more theoretically sound ECRA algorithm yield good results? In figure 11 we show that actually the ECRA BFM algorithm did yield good results, and identifies the preclusters properly, even finding the proton cluster that broke down.

With these three algorithms in mind, we build three different cluster recognition tools: MSTE, ECRA and MSTpC. MSTE and ECRA are the regular algorithms, while MSTpC is the proton MST algorithm with the cloud of neutrons that are near each MST cluster. In figure 12 we show the evolution of the size of the largest fragment for the early stages of the evolution for the three clusters: MSTE, MSTpC and ECRA. The figure shows that the ECRA fragment remains relatively stable and stabilizes quickly, while the other two algorithms yield fragments that are always larger and stabilize more slowly. It is also interesting to note that the MSTpC fragment starts with about 100 neutrons more than the corresponding ECRA fragment, which means that the ECRA fragment is in a very neutron rich environment. This kind of situation makes the r-process more likely to happen.

On the other hand, in the expansion of the lasagna structure, none of the algorithms for fragment recognition identifies clusters very early on the evolution (see figure 13). At this stage, there is a very large fragment, which is actually infinite. Nevertheless, the ECRA analysis shows that this fragment breaks down early into many different fragment and, as a result, the mass of the largest fragment decreases drastically with time. It is interesting to notice that, unlike the pregnocchi, in this case the MSTpC algorithm is the one that takes the longest to identify that the infinite cluster breaks down. This shows us that the ECRA algorithm is also more verstatile to study the early fragment formation. It is also of interest (not shown) that the proton fraction $x$ of these fragments is relatively stable for the ECRA algorithm, while the other two yield a proton fractio that decreases monotonically with time.

VI Discussion and Concluding Remarks

We studied, with molecular dynamics, structural properties of the crust of neutron star through three different potentials. These potentials involve a nuclear term tailored to reproduce binding energies and compressibilities of nuclear matter plus screened Coulomb interaction. To analyze the structures formed, we used four different cluster recognition algorithms: MST, MSTE, MSTpC and ECRA-BFM. With these algorithms we found that of the three potentials, two of them (New Medium and SSP) developed a newly found structure for low proton fractions that we called pregnocchi. This structure consists of proton aggregates formed by the mediation of the attractive $V_{np}$ term of the potential that withstood the expansion.

We also analyzed the expansion of the infinite neutron rich matter described in terms of the little big bang model. We showed that in general the proper identification of the structure is highly dependent on the algorithm chosen, being ECRA and MSTpC the most suitable to find the structures and ECRA the most stable one. This approach, combined with different cluster algorithms, allowed us to identify the dynamics of the fragment formation. The asymptotic state showed a high dependence on the rate of expansion, both in the mass histogram and the spatial distribution of the fragments: for fast enough rates, the expansion was similar to a mechanical fracture, where the spatial distribution was heaviliy correlated with the original. However, for slower rates, the expansion was a thermal expansion in which the asymptotic state was relatively homogeneous. The clusters formed in the slower expansion were much larger than those formed in the fast expansion. A thorough analysis of the clusters formation dynamics showed that they were formed early in the expansion. In particular, the novel structure that we have called pregnocchi is quite relevant, because according to ECRA analysis these preexistent aggregates evolve in time embedded in a neutron cloud, giving rise to configurations in which r-proccess might set in.

Appendix A On the stability of MSTE clusters

A simple example can be studied to see whether MSTE clusters are always stable. Consider an interaction

[TABLE]

Now we study a set of particles of mass $m$ with positions $r_{i}=i\,a$ (with $i\in\mathcal{Z}$) so that every particle is at a distance $a$ from its nearest neighbors. If the velocity is $v_{i}=i\,v$, each particle will be energetically bound with its neighbors if $v\leq\sqrt{2\,V_{0}/m}$. For $2n+1$ particles, with $-n\leq i\leq n$, the kinetic energy of the system will be

[TABLE]

The potential energy, however, is

[TABLE]

It is clear then that for large $n$, no matter the value of $v$, the system will be unstable even though MSTE algorithm recognizes it as a single cluster.