Three dimensional structure of low-density nuclear matter

We numerically explore the pasta structures and properties of low-density nuclear matter without any assumption on the geometry. We observe conventional pasta structures, while a mixture of the pasta structures appears as a metastable state at some transient densities. We also discuss the lattice structure of droplets.


Introduction
In low-density nuclear matter, which is realized in the supernovae core or in the crust of neutron stars, inhomogeneous structures of nuclear matter are expected [1]. With increase of density, which ranges from well below to the normal nuclear density, the shape of constituent nuclei is expected to change from spherical droplet to cylindrical rod, slab, cylindrical tube, spherical bubble, and to uniform. These shapes are figuratively called as "nuclear pasta".
The density-dependence of the species and the sizes of pasta are determined by minimizing the total energy density, i.e. the sum of the bulk, the surface, and the Coulomb energy densities. We can roughly assume that nuclear matter at sub-saturation density consists of a dilute gas phase and a dense liquid phase in chemical equilibrium, which determines particle densities in both phases. Thus the bulk energy density is independent of the structure, once the volume fraction is given, but the shape and the size of the structure are determined by the balance between the surface and the Coulomb energy densities.
Among the early studies of nuclear pasta, geometrical symmetry of the structure was very often assumed, i.e. the Wigner-Seitz (WS) cell approximation was employed: The whole space is divided into equivalent cells with charge neutrality, and a geometrical symmetry is assumed with a given dimensionality. Then the shapes of the cell becomes sphere in three dimension (3D), cylinder in 2D, and plate in 1D cases. Reflective boundary condition is imposed to the density distribution of each particle. There is no interaction between cells due to the charge neutrality. So all the physical quantities of matter are represented by those of a single WS cell. Furthermore, the calculation is only in one-dimension due to the symmetry, which drastically reduces the computation.
Within this approximation and assuming uniform density distributions in each phase, analytic expressions of the surface and the Coulomb energies can be derived for a given dimensionality of the symmetry. In 1983 Ravenhall et al. have shown [2] that the configuration to give the lowest energy changes from spherical nuclear droplet (3N) to cylindrical nuclear rod (2N), nuclear slab (1NB), cylindrical bubble (2B), spherical bubble (3B) and to uniform (U) with increase of the density. In 1984 Hashimoto et al. have included the Coulomb interaction among cells [3] and have got essentially the same results. In this work the Coulomb energy was evaluated without the WS approximation, but only simple shapes of nuclei were assumed.
On the other hand, Williams and Koonin have calculated nuclear pasta using the Thomas-Fermi model for a system in a cubic cell with periodic boundary conditions [4]. In this calculation no assumption was made for the structure of matter and they got essentially the same results as the studies with the WS cell approximation. There is also a recent calculation with the Hartree-Fock theory in a cubic cell with periodic boundary conditions, giving almost the same results again [5]. Though the above calculations did not assume any particular structure, the size of the cell was rather small so as to include only one period of the pasta structures. The usage of small cell brings about an implicit but strong restriction to the matter structure, suppressing the appearance of complex structures.
In this article we perform three-dimensional calculations in periodic cubic cells with sufficiently large sizes. We will discuss how the pasta structures appear in the ground state of matter, showing their crystalline structures. We also show some metastable states of matter which may be realized at finite temperatures.
Equations of motion for fermions simply yield the standard relations between the densities and chemical potentials, where m * N (r) = m N − g σN σ(r) is an effective mass. We assume here the global charge neutrality. In case of cold catalyzed matter, the relation µ n = µ p + µ e is further imposed.
To numerically simulate infinite matter, we use a cubic cell with periodic boundary conditions. We divide the cell into three-dimensional grid points. The best cell size should be as large as to include several periods of the pasta structures, and the best grid width should be as small as to attain smooth meson and fermion fields. Due to the limitations in the memory and the CPU time, we use the cell size of ∼ 60 fm and the grid width |dr| ∼ 0.8 fm. If only one or two periods of structure appear in the calculation cell, its shape may be affected by the boundary condition and the appearance of some structures is implicitly suppressed.
Giving average densities of baryons (ρ B ) and electrons, initial density distributions are randomly prepared. Then proper density distributions and meson fields are searched for. To obtain the density distributions of baryons and electrons we introduce the local chemical potentials µ a (r) (a = p, n, e). The equilibrium state is determined so that the chemical potentials are independent of the position. An exception is the region with an empty particle density, where the chemical potential of that particle becomes higher. We repeat the following procedures to attain uniformity of the chemical potentials. A chemical potential µ i (r) of a baryon i = p, n on a grid point r is compared with those on the six neighboring grids r ′ = r + dr, (dr = ±dx, ±dy, ±dz). If the chemical potential of the point under consideration is larger than that of another µ i (r) > µ i (r ′ ), some part of the density will be transferred to the other grid point. This adjustment of the density is done simultaneously on all the grid points. In addition to the above process, we adjust the particle densities between distant grid points chosen randomly so as to avoid making separate droplets with different µ i .
The meson fields and the Coulomb potential are obtained by solving Eqs. (1)-(4) using the baryon density distributions ρ i (r) (i = p, n) and the charge density distribution ρ p (r) + ρ e (r). These equations are solved numerically by a conjugate gradient method. The electron density ρ e (r) is directly calculated from the Coulomb potential V Coul (r) and the electron chemical potential µ e . The global charge neutrality is then attained by adjusting µ e . Above processes are repeated many times until we get convergence. We used PRIMERGY BX900 of JAEA massively parallel computing system. To finish a calculation of a typical case, about 10 CPU days is needed.

Result
We present here our first results with fixed proton fraction Y p for Y p = 0.1, 0.3 and 0.5, leaving catalyzed matter in another paper. Shown in Fig. 1 are the proton density distributions in cold symmetric matter (proton fraction Y p = 0.5). We can see that the typical pasta phases with rod, slab, tube, and bubble, in addition to the spherical nuclei (droplet), are reproduced by our calculation in which no assumption on the structures was used.  Any intermediate structures does not appear as a ground state at any density. In a droplet, we have seen that the proton density is highest near the surface due to the Coulomb repulsion, while the neutron density distribution in the droplet is flat. Note that baryon density outside the droplets is zero for Y p = 0.3 and 0.5. The electron density is finite over all space but slightly localized around the droplets. We can see this behavior of fermions for Y p = 0.5 and ρ B = 0.016fm −3 in Fig.  2, where plotted are the densities of proton, neutron and electron along a line which crosses through the droplets.
We show the density dependence of the energy, the total pressure, and the baryon partial pressure in Fig. 3. The density dependence of these pressures is qualitatively the same as the one with the WS cell approximation. The difference between our results and those with the WS cell approximation, in which the same RMF framework is used, is the density region of each pasta structure; density region of the rod is wider and the tube narrower in our calculation. Figure 4 shows the radius of droplets R d , the lattice constant R cell , and the volume fraction u of droplets. Here, R d and R cell are defined as follows, and u = (R d /R cell ) 3 , where V denotes the cell volume, N d the number of droplets in the cell, and the bracket ... the average over the cell volume. Comparing the present results and those with the WS cell approximation, the volume fraction shows the same behavior, but the lattice constant and the radius of the droplets are different. When a spherical nucleus becomes large, it becomes unstable against fission. It is expressed as the Bohr-Wheeler condition [7], Coul is the Coulomb energy of a nucleus. In the WS cell approximation, the Coulomb energy in a cell is expressed using the Coulomb energy of a nucleus as E Coul ≈ E (0) Coul 1 − 3u 1/3 /2 . On the other hand, the condition of optimum nuclear size gives E surf = 2E Coul . From these equations, the appearance of non-spherical nucleus in nuclear matter due to the fission instability has been expected for the volume fraction u > 1/8. However, in our calculation, the structural change from droplet to rod occurs around a volume fraction u = 0.2 (see the value at ρ B ≈ 0.02 fm −3 in Fig. 4). The relation between the Coulomb energies of a cell and that of a nucleus has been derived by using a uniform background electrons and uniform baryon density in a nucleus. The effect of the screening by charged particles, which is naturally included in both the WS cell calculation in Ref. [6] and our present calculation, may be one of the origins of this difference. Also the difference of the droplet surface may be the origin, since they used the compressible liquid drop model with a sharp surface, while our droplets have a diffuse surface.
There are some differences between our calculation and the one with the WS cell approximation. Calculating in a large cell which has several periods of pasta   Figure 4: Density dependence of the radius, the lattice constant and the volume fraction. Red lines are our results and green lines are those with the WS cell approximation One point unlike the conventional results emerged in the crystalline structure of droplets. In our calculation, it emerges as a face-centered cubic (fcc) lattice, while it has been regarded to take a body-centered cubic (bcc) lattice in the previous studies [8,9]. Crystalline structures in the bcc and fcc lattices give rise to a subtle difference of the Coulomb energy, which amounts to about 0.2-0.8 MeV. The ratio of the Coulomb energy for total energy difference is about 20%. In other words, most of energy differences come from the bulk energy of the nuclear matter.  Table 1 shows the baryon-density dependence of radii of droplets in the cases of the fcc and bcc lattices. The radii of droplets are different even if their baryon densities are the same. In Refs. [8,9], it is reported that the bcc crystalline structure of droplets is realized in the ground state at low densities, while the fcc crystalline structure appears in our calculation. This difference may partially come from that they compare the bcc and fcc crystalline structures using droplets of the same size, while in our case they have different droplet sizes. In the QMD calculations [10,11] that precede the present calculations without assuming structure not only for static but also for dynamical aspects of nuclear matter, droplets form the bcc crystalline structure. This difference might come from the treatment of electrons: In the QMD calculation, uniform electron distribution has been assumed, while in our calculation, inhomogeneous electron distribution is attained consistently. We have performed another calculation with uniformly distributed electrons. However, the result did not change and the fcc structure is more favored. Thus it is confirmed that the charge screening by electrons does not very much affect the crystalline structure [6]. This difference of the crystalline structure between the QMD calculation and the present calculation remains as a future problem.
In Figs. 5 and 6, we show the density dependence of the energy, the total pressure and the baryon partial pressure for Y p = 0.3 and 0.1, which are roughly relevant to the supernova matter and the neutron star crust. We obtain the typical pasta structure as ground states for any proton fraction above 0.1. Also in the cases of Y p =0.3 and 0.1, the fcc lattice of droplets is energetically more favorable than the bcc. For rod and tube, the simple lattice is more favorable than the honeycomb lattice. In Fig. 7 we plot the density profile of proton and neutron for Y p = 0.1 with baryon density 0.016fm −3 along a line which crosses through the droplets. The neutron density is finite at any point: the space is filled by dripped   neutrons, which is in contrast to the case with higher Y p where the neutron density is zero outside the nucleus.   On the way of searching for the ground-state structures, we sometimes observe exotic structures which are energetically metastable. Figure 8 (a) shows a mixed structure of droplet and rod at ρ B = 0.022 fm −3 . Similarly, (b) is a mixed structure of slab and tube at 0.068 fm −3 . These structures appear around densities where the structures change. If these structures had appeared as ground states, transition from droplet to rod, and from slab to tube would happen more smoothly by way of those intermediate structures. In fact, the QMD calculations have reported some intermediate structures of droplet and rod, slab and tube as ground states [10].
We have observed more exotic structures: Panel (c) at 0.018 fm −3 looks like dumbbells which have been reported to appear in the dynamical compression of matter [12]. Finally, panel (d) at 0.044 fm −3 is a diamond structure, which resembles double-diamond structure studied with a compressible liquid drop model [9,13]. These complex structures are difficult to be the ground states since they have larger surface area than simple pasta structures. At finite temperatures, however, these structures might contribute to the Boltzmann ensemble.

Summary
We have numerically explored inhomogeneous structures and properties of low-density nuclear matter using the RMF model and the Thomas-Fermi approximation. Without any assumption on the geometry, we have carried out fully three-dimensional calculations on large cubic cells with periodic boundary conditions. With increase of density, which ranges from well below to the normal nuclear density, we have observed that the ground state of matter shows the typical pasta structures. More complex structures like "diamond" structure, "dumbbell" structure, and mixtures of two types of pasta structures appear as metastable states at some transient densities. As for the crystalline structures for the droplet and the bubble structures, fcc lattice has been more suitable than bcc lattice, which is different from the previous studies. For neutron star crust and supernovae matter, we should explore these structures and properties for β-equilibrium and extend to finite temperature. At finite temperatures, complex structures might contribute to the Boltzmann ensemble. We can apply these contributions to the mechanism of glitch and cooling process of neutron stars and thermal and mechanical properties of supernovae. To explore β-equilibrium matter, we need more calculation space. But there are some problems in calculation performance. As mentioned before, to finish a calculation of a typical case, many CPU time is needed. Most of the CPU time is consumed by the part to get the uniformity of chemical potentials and its iteration. To reduce the CPU time, we must improve the method and calculation code of this part. More realistic calculation can be done by, for example, including gradient terms of the density to improve the description of the surface properties [14].