Nuclear structure calculations for neutron-star crusts
The goal of this paper is to investigate properties of clusterized nuclear matter which is believed to be present in crusts of neutron stars at subnuclear densities. It is assumed that the whole system can be represented by the set of Wigner-Seitz cells, each containing a nucleus and an electron background under the condition of electroneutrality. The nuclear structure calculations are performed within the relativistic mean-field model with the NL3 parametrization. The first set of calculations is performed assuming the constant electron background. The evolution of neutron and proton density distributions was systematically studied along isotopic chains until very neutron-rich system beyond the neutron dripline. Then we have replaced the uniform electron background with the realistic electron distributions, obtained within the Thomas-Fermi approximation in a self-consistent way with the proton distributions. Finally, we have investigated the evolution of the -stability valley as well as neutron and proton driplines with the electron density.
description=average electron density
description=average baryon density
description=average baryon density
description=electron chemical potential
description=difference of electron chemical potentials
description=electron chemical potential
description=electron chemical potential of uniform electron distribution
description=neutron chemical potential
description=proton chemical potential
description=normal nuclear density
description=total charge density
description=kinetic energy of uniform electron distribution
description=separation energy of protons
description=total nuclear energy
description=total nuclear energy of uniform electron distribution
description=Wigner-Seitz cell energy
description=difference of Wigner-Seitz cell energy
description=kinetic energy density
description=Wigner-Seitz cell filling fraction
description=number of baryons
description=number of neutrons
description=number of protons
description=Wigner-Seitz cell volume
\newglossaryentryrad_wsc name=, description=Wigner-Seitz cell radius
\publishers1) Frankfurt Institute for Advanced Studies, J.W. Goethe
University, Ruth-Moufang Strasse 1,
60438 Frankfurt am Main, Germany
2) National Research Centre Kurchatov Institute, 123182 Moscow, Russia
The properties of clusterized nuclear matter at subnuclear densities and moderate temperatures are under active investigation for several decades, for a review see refs. [1, 2, 3]. After the seminal papers [4, 5], this topic became a very active research field in recent years. Many microscopic and macroscopic models have been used, see e.g. refs [6, 7, 8, 9]. The superfluid properties of such inhomogeneous matter have been studied in refs [10, 11, 12]. Often the calculations are done within the Wigner-Seitz approximation, when the system is split in electrically-neutral cells containing a nucleus and electrons.
In our previous paper , we have studied the influence of a uniform electron background on nuclear ground states and decay modes. The nuclear structure calculations were performed within the Relativistic Mean Field (RMF) model with parameterization NL3 . It was found that for large electron Fermi momenta, decay and spontaneous fission are significantly altered, and, depending on the electron density, lead to stabilizing or destabilizing effects. With growing electron density, the -stability valley moves to more neutron-rich systems, and the proton drip line moves to more proton-rich nuclei. These changes happen because the single-particle potential for protons gets deeper as electrons are added. Since neutrons are electrically neutral, in the first approximation the neutron drip line is unaffected by the presence of electrons. It was found that for electron Fermi momenta of , the -stability line reaches the neutron drip line, and free neutrons occur. Hence, for larger electron densities, these quasi-free neutrons need to be taken into account.
In the present paper, we extend this model in two directions. We include the dripping-out neutrons which populate unbound states. As before, we adopt the Wigner-Seitz approximation representing the whole system as a collection of non-interacting electrically-neutral cells characterized by the nuclear charge , baryon number , the ratio and average baryon density , where and is the WS cell radius. We employ boundary conditions for the single-particle states which are physically motivated and which provide a non-zero density of the neutron gas at the cell boundary.
In  we have shown that using realistic electron distributions obtained by solving the Poisson equation leads to a significant reduction of the electrostatic energy and electron chemical potential. In the present work, we have calculated the electron density distributions consistent with the proton distributions obtained self-consistently within the RMF model. We calculate the whole nuclear chart up to very heavy nuclei which can be formed in neutron star crusts. Both nuclear matter and the electron distributions in a WS cell are calculated using an iterative algorithm. We demonstrate that the realistic electrostatic potentials induce a shift of the -stability valley and the proton drip line to more proton-rich systems thus utilizing the stabilizing effect of the electrons. The differences in the electron chemical potentials are evaluated for different values of and the baryon density.
2 The Framework
2.1 The RMF model
The nuclear structure calculation below are carried out within the same RMF-NL3 model as in ref. , but now the unbound neutrons are allowed too. The corresponding Lagrangian has the form
where denotes the nucleons (neutron or proton) and are the occupation probabilities of the nucleon states. They take into account the BCS pairing with a density-indepentent -force . This approach allows us to treat both bound and unbound nucleons on the same level. In this paper we assume zero temperature.
2.2 The Wigner-Seitz approximation
The calculations below are performed within the Wigner-Seitz (WS) approximation by dividing the system into WS cells, each containing only one nucleus and the number of electrons equal to the nuclear charge . The introduction of WS cells is aimed as a convenient approximation of an ensemble of nuclei surrounded by a more or less uniform background of electrons and neutrons. Its construction should ensure that the physics contained in an individual cell is to large degree independent of the surrounding matter. This requires, for example, charge neutrality of the cell. However, there could be also other requirements which minimize the electrostatic interactions with neighbouring cells (see below). The accuracy of the WS approximation for the inner neutron-star crust has been discussed, e. g. in Refs. [17, 18]. We use the units with .
2.3 Boundary conditions for unbound neutrons
In our calculations we allow the possibility that some neutrons are not bound in a nucleus, but occupy the whole available volume and reach the boundary of the WS cell, which in our scheme coincides with the edge of the numerical grid. For this reason, boundary conditions for the neutron states at the cell edge need to be specified. The choice of these boundary conditions should be guided by physical considerations. Physically, the unbound neutrons are expected to form a more or less homogeneous gas, in which the nuclei are embedded. Thus, the neutron density at the edge of the cell should have a value which coincides to a large degree with the average density of unbound neutrons. Futhermore, the slope of the neutron density should be zero in order to have a continuous gas in the region between the cells and to allow a continuous connection of physical quantities across cell boundaries.
In the literature, different possibilities have been discussed. As demonstrated in , reasonable properties of the neutron states can be achieved if the boundary conditions are chosen differently for states with odd and even angular momenta. For spherical cells, such a a sorting of states is possible since the angular momentum is a conserved quantity in a spherically symmetric system. One possible choice is to set
This way, approximately half of the states assumes a zero value at the boundary, while the other half assumes a zero slope. There is no physical reason why the condition above should be chosen this way and not the other way round, i.e., that states with odd angular momentum are required to have zero slope and states with even angular momentum are required to have zero value. Baldo et al  have investigated the effect of a different choice on physical oberservables and found a non-neglegible effect.
As has been shown already in , these boundary conditions are often not sufficient to prevent a building up of density at the cell edge during the numerical iterations. One possible reason for such a buildup during the iterative procedure is given by the fact that – as many neutron states are involved – many states have large angular momenta and are thus located at rather large radii. By this reason, the density accumulates near the cell edge, giving rise to an attractive mean field potential for other neutron states. Therefore, during the iterations, more and more neutron states are attracted to the border of the cell.
In , this problem has occured in the framework of non-relativistic mean-field calculations. It was cured by an averaging procedure; i.e. both the particle density and the kinetic energy density were averaged over a radial region in the cell where the neutron gas density was approximatey constant. The densities in a region close to the grid egde with were then substituted with these values. In each iteration, the density close to the cell edge was brought back to the value of the average neutron gas in the cell
We checked this averaging procedure in our calculations and made several test runs. While the procedure worked in principle, it lead to some undesired effects. For some nuclei and cell radii, this procedure led to a well-behaved, smooth and almost constant neutron gas within the cell. For certain cases, however, the averaging near the cell edge cured the buildup of density in each iteration but did not lead to a reduction of the density at the cell edge during the iterative procedure. In each new iteration density was still summing up at the cell edge. If such a condition remained until the convergence of the solution, the net effect resulted in a removal of neutrons from the cell, since the large rise of the neutron density at the cell edge was replaced by average values, and no measure was taken to conserve neutron number. In extreme cases, a reduction of the neutron number by up to was found. These artifacts showed up in observables such as the binding energy of the system. Since the particle number in spherical calculations scales with , even small changes of the density at the boundary can lead to a large error.
To conserve the particle number and get rid of these side effects, a different method was designed and implemented in our present calculations. The behavior of the neutron gas near the cell edge has been constrained by the physical condition that the neutron density should have a zero slope at the cell edge. Such constraints in self-consistent calculations can be used in all kinds of situations where some nuclear properties should have prescribed values. For example, the fission barrier of a heavy or superheavy nucleus can be computed by putting a constraint on the quadrupole moment of the nucleus and by varying its value. In our case, the constraint was chosen such that the neutron gas density was required to have zero slope at the boundary. It was implemented by simply requiring that the difference of the values of the neutron density at the two outer-most grid points should be zero, i.e.
where is the radial spacing of the numerical grid in coordinate space. In principle, more constraints could be added regarding other pairs of grid points near the border of the grid, but with respect to the desired speed and convergence of the calculations this most simple version proved to be most adequate.
In fact, this constraint can be applied not only to the neutron density but to the total baryon density as well. In most cases, the protons are located at the cell center and form bound states. Then the proton wave functions are approximately zero at the cell edge, and the proton density does not contribute to the total baryon density there. In cases when the protons dissolve and the proton states reach out to the cell edge, such a requirement is fully justified. The use of this constraint in our calculations led to satisfactory results in most cases. The particle number was always exactly conserved. Only for very large neutron numbers and/or comparatively small cell radii or large average values of the neutron density, a small density bump close to the cell edge appeared. The constraint was employed during the whole iterative procedure or switched on at a later stage, depending on the system size and the convergence speed.
2.4 Initialization of states
For the considered extremely neutron-rich nuclear systems, the initialization proved to be quite important. During a self-consistent iteration of a regular bound nucleus, more states than nucleons in the system are usually initialized, and the lowest states are occupied with nucleons. During the iterative procedure, however, non-occupied states can become energetically more favorable than the occupied states, in which case these states will be occupied with nucleons. For systems with a moderate number of states, this scheme works extremely well and provides the configuration with the lowest ground-state energy.
For systems with many neutrons, a straight-forward initialization of all states from scratch can be problematic. Firstly, the states are usually initialized in a harmonic oscillator well, which has only bound states. Thus, these states are quite different from the unbound states which they will become during the iterative procedure. Therefore, rapid changes of these states will occur during the first few iterations. Secondly, for the ground-state configuration, many states should be bound within the nuclear potential.
For these reasons, for the calculation of a nuclear system with protons and neutrons, where , the following iterative procedure was used: in a first step, a nucleus with protons and neutrons was calculated, where corresponds to a neutron number for which the nucleus is close to the line of stability. After the solution has converged, the proton and neutron states (both occupied and unoccupied) are stored. In a second step, these states are used for the initialization of a nucleus with protons, but neutrons, where . Relations or proved to work well. After such a chain of calculations, a well-behaved solution for the extremely neutron-rich system with protons and neutrons could be obtained.
2.5 Simplified treatment of electrons
Following our previous work  for rough estimates, we employ in this section a simplified model assuming a uniform electron density. Generally, the chemical potential of the electrons can be represented as
The kinetic part is simply given by
where the electron Fermi momentum does not depend on . The second contribution to the chemical potential stems from the electrostatic potential in which they move. This potential is produced by the total charge density , which is given by both electrons and protons
The electrostatic potential is thus space-dependent due to the spatial dependence of the proton density. Thus, for given constant and space-dependent , would have to be space-dependent, too, which is a contradiction to the assumption of a constant electron background. This inconsistency can be cured by some extent by introducing a constant potential , defined as the average electric potential felt by an electron.
Since is defined as
the average electric potential can be calculated as
where is the total number of electrons (and due to charge neutrality) equal to the number of protons in the system. Thus, in the calculations presented in this section we use the expression
The condition of -stability is given by
For zero temperature, the Fermi energies of protons and neutrons correspond to their respective chemical potentials. The Fermi energies are obtained from the iterative calculation. When pairing is employed, the Fermi energy is not uniquely determined but lies somewhere between two states such, that the average nucleon number is fulfilled. Due to these uncertainties, we formulate this condition as
3.1 Evolution of density profiles
First we investigate the structural evolution of shapes along isotopic chains. It is assumed that the nucleus is located at the center of the spherical Wigner-Seitz cell of radius . We start with a -stable isotope and add neutrons. In Fig. 1 left and right shows, respectively, the neutron and proton densities for tin isotopes. The drip-line isotope is marked by crosses. These isotopes correspond to different total baryon numbers within the cell and thus to different baryon densities.
Up to the drip-line isotope with , the neutron density shifts to larger radii, but all neutrons are bound. Beyond the drip line, neutrons leak out from the nuclear potential and start to form a constant neutron gas filling the cell approximately uniformly. For , the neutron density at the cell edge is approximately . The corresponding average total baryon density in the cell is . These numbers are practically equal because i) and ii) the dripped neutrons are distributed more or less uniformly in a big volume. Also note that the central neutron density remains almost constant for all isotopes, , the value corresponding to the saturated nuclear matter.
For all neutron numbers, the protons remain localized within the central part of the cell. With increasing neutron number, however, their distribution broadens. In contrast to the neutron distributions, the central density of the protons changes with increasing neutron number. While it is for the beta-stable isotope, the central proton density drops below for the heaviest isotope.
In our opinion, this behavior can be explained by the symmetry energy contribution. Indeed, when the dripped neutrons build up a significant density outside the nucleus, it is energetically favorable to stretch the proton distribution in order to reduce the neutron-proton asymmetry in the transition region. This stretching effect is rather robust. It is clearly seen in all three considered cases. According to our knowledge, this effect was never discussed previously. One should also bear in mind that the increased asymmetry in the central region is less important in the energy balance because of the weighting of different spherical shells.
Analogous plots for lead isotopes () and isotopes of the superheavy element are shown in Figs. 2 and 3, respectively. The neutron distributions in both cases are rather similar to the one for Sn isotopes. But the proton distributions show large differencies. In the case of , the protons have a small bump at the center. But in the case of the superheavy system , protons have a pronounced central deep, which goes down to for . According to our knowledge, such hollow shapes have not been discussed previousely. In our opinion, this is a consequence of the strong Coulomb repulsion of protons in such a heavy system.
3.1.2 Constant electron to baryon fraction
In this subsection we present calculations for sequences of nuclei with constant proton-to-baryon ratio
which corresponds also to the fixed electron fraction in the charge-neutral Wigner-Seitz cell. Such constraint is often used to specify physical conditions in stellar matter.
In Fig. 4 we show the baryon density profiles for and an average baryon density in the cell of , which is about 75 times lower than the normal nuclear density. For (the bottom panel in the figure), the corresponding neutron number is . The top panel corresponds to proton number and neutron number . For each panel, the dashed line indicates zero density. One can see that the profiles evolve slowly from the Woods-Saxon type of distributions at to more hollow shapes at larger .
This trend becomes even more pronounced for higher baryon density and smaller , as shown in Fig. 5. In this case the hollow shape gradually evolves to the extreme configuration where nucleons are located on a spherical shell with a width of about . As one can see, the central density in this case is not much different from the density near the WS cell boundary. Again, the dashed line indicates the zero density level. In Fig. 5, the proton number and neutron number are chosen for the bottom panel and proton number and neutron number are fixed in the top one. One should bear in mind that these nuclei are generally not in -equilibrium.
3.2 Single-particle states
3.2.1 -stable isotopes
Let us consider now the -equilibrium conditions in WS cells which contain unbound neutrons and an uniform electron background. At fixed average baryon density , such systems can be specified by the total charge (isotopes) or baryon number (isobars). Then, the volume of the cell is obtained self-consistently by the iterative procedure. In Fig. 6, the occupied single particle states for protons (left) and neutrons (right) are shown within the respective single-particle potential for the system with , which is the -stable tin isotope for the average baryon density . The arrow shows the electron chemical potential, which is added to the proton chemical potentialto reach the upper occupied level of neutrons, see eq. (11).
The occupied proton states are well bound in a potential well that resembles the proton potential for ordinary nuclei, with the exception that it does not go to zero as increases but reaches a constant (negative) value. This is due to the constant neutron and electron backgrounds which generate an attractive potential. The lowest neutron levels are well bound in the potential well. At the energy level of about , the neutron potential becomes flat up to the cell boundary. At this point, a drastic increase of the level density is observed.
While for ordinary nuclei the proton potential is less deep compared to the neutron potential, we encounter an opposite situation here. Since we require charge neutrality of the Wigner-Seitz cell, the cell is filled with electrons which interact attractively with the protons, and do not directly interact with the neutrons. There is an indirect effect due to the self-consistency of nuclear interactions, but it is rather small. Hence, due to the attractive electron background, the repulsion of the protons is reduced compared to the situation in vacuum.
The -stable tin isotope that corresponds to the baryon density is shown in Fig. 7. Since the baryon density is larger than in the previous case, the radius and correspondingly the volume of the WS cell is smaller. Hence the electron chemical potential is larger, and in order to reach -equilibirum, more neutrons should be added to the cell. Consequently, neutrons in this nucleus occupy higher-lying states.
Due to the selfconsistency of the calculation, especially due to the interaction of protons with neutrons and electrons, the proton single-particle states in the WS cell with , have shifted down, but the neutron states have shifted up. The general picture is not very different for a heavier system with , shown in Fig. 8.
3.2.2 -stable isobars
It is interesting to see how the nuclear structure changes with increasing baryon density for isobars, i.e. at fixed and different . Two isobars of are presented in Figs. 9 and 10 for baryon densities and , respectively. Since the electron chemical potential rises with the density, more protons convert to neutrons at higher baryon densities to maintain the -equilibrium. As one can see, indeed at the larger baryon density shown in Fig. 9, the system has larger and smaller compared to the one shown in Fig. 10.
4 Nuclear structure calculations with realistic electron background
4.1 Electron density distributions
In this section, we extend the model by replacing the constant electron density with a realistic, nonuniform density distribution obtained self-consistently by using the relativistic generalization of the Thomas-Fermi method. In the quasiclassical approximation one can write the energy of an electron of momentum at point as
All states with are occupied, where is the electron chemical potential independent of . The local Fermi momentum is determined by the condition
Summing all the electron states from up to the local Fermi momentum we get the local electron density
where is the spin degeneracy factor. The requirement of electroneutrality of the cell implies that the total electron number is equal to the total proton number in the cell volume ,
which determines the electron chemical potential .
The proton density is obtained from the RMF calculation, so the electrostatic potential can be determined self-consistently from the Poisson equation
where is the total charge density defined in eq. (7) which itself depends on via eq. (17). This equation is solved with the boundary conditions for requiring regularity at the origin and vanishing derivative at the cell boundary. The nonuniform electron density distribution calculated in this way determines the electrostatic potential for the next iteration of the RMF calculation. All particle densities, including electrons, are obtained self-consistently. In order to use the same notation as in the previous sections, we define (without argument ) to characterize the average electron density in the WS cell,
4.2 Numerical implementation of self-sonsistent approach
Now we are going to study the properties of nuclei embedded in the realistic electron background, in particular, their binding energies as well as neutron and proton drip lines. To fulfill this task one should evaluate a broad ensemble of particle-stable nuclei, i. e. the whole nuclear chart. Certainly, this task presents a serious computational challenge. Below we limit ourselves to the nuclei with even proton numbers from 8 to 240 and neutron numbers from 6 to 360. Then the total number of WZ cells which should be coinsidered amounts to about 60000. Each individual cell calculation takes a CPU time between minutes and hours, depending on the electron Fermi momentum and the nucleon number. For these extensive computations we have used the CPU cluster at LOEWE Center for Super-Computing (CSC) at Goethe University. We carefully selected appropriate configurations for the SLURM workload manager used at this cluster and wrote special scripts to assign the amount of cells to allow parallel computation.
Auxiliary scripts created the input files for each WS cell system calculation and assigned the calculations to the corresponding computer nodes. The iteration scheme for each system started with the calculation of the baryon densities within the RMF description, where the solution of the Dirac equation gave the single-particle wave functions. The many-particle wave function was then obtained by calculating the Slater determinant. We used cubic B-Splines to interpolate between the numeric grid of the nucleons and the electrons. The Thomas-Fermi routine calculated the self-consistent electrons by iterating between the solution of the one-dimensional Poisson eq. (19) and the electron density distribution given by eq. (17). The numerical scheme of the Thomas-Fermi routine is described in .
4.3 Modification of nuclear chart under stellar conditions
In this subsection, we present results of our nuclear structure calculations for characteristic conditions in outer crusts of neutron stars. They can be conveniently specified by the average electron density or the Fermi momentum as defined in eq. (20). We believe that rather heavy nuclei can exist in the crusts of neutron stars or supernova interiors at subnuclear densities and low temperatures. Numerous calculations within the statistical approach predict broad distributions of nuclei up to mass numbers 400-500, see e.g. recent paper  and references therein. As was demonstrated in our previous work, the spontaneous fission of such nuclei is suppressed because of the screening effect of electrons.
Our particular interest was to compare results obtained in  within a simplified model of uniform electron background with more realistic calculations presented in this paper. To better illustrate the difference, we calculate the energy per nucleon defined as
where is the total energy of the nucleus (with respect to the rest mass mass of nucleons), and is the change of the electron kinetic energy between realistic (r) and constant (c) density distributions in each cell. This change is significant for heavy nuclei, as we have shown in .
In Fig. 11, the results are presented for the electron Fermi momentum of . One can see that the energy gain per nucleon ranges from for light to for heavy nuclei. Beyond the proton dripline, in addition to nuclei and electrons, the cells will contain also free neutrons. It is interesting to note, that the valley of the most bound nuclei in general does not coincide with the beta stability line, defined by the condition on the electron chemical potential .
In Fig. 12 we present the proton drip lines for a wider range of electron densities from to . They are compared with the previous calculations  for constant electron densities. The drip lines are defined in terms of the separation energies which is approximated by a half of the two-proton drip line, . By definition, the flip of the sign marks the position of the dripline.
As expected, the effect of realistic electron distributions is strongest for most heavy nuclei and highest Fermi momentum of electrons. For instance, for , we have found for a shift of the proton dripline from to . At smaller , the shift varies from 1 to 2 units of charge. Still, even for , the shift can reach a few charge units, e.g. for from to , as shown in Fig. 13. Similar shifts are also predicted for the charge and mass numbers of most bound nuclei.
As demonstrated in ref. , using realistic electron distributions is very important for calculating the electron chemical potential . The correction depends on the average baryon density as well as the electron-to-baryon fraction . Figure 14 shows the shift of the electron chemical potential for five different proton numbers (\glsnum_prn , , , and ) with two isotopes each, calculated for varying electron Fermi momenta up to .
Following ref. , we define the cell filling fraction as
where is the mean baryon density und is the normal nuclear density as in . In order to compare our results with our previous calculations in ref. , we show the values of the cell filling fraction (in percent) for each isotope curve on the right side of Fig. 14. The numbers are given by the right side of eq. (22) for and for (as used in refs.  and ). Obviously, the electron-to-baryon ratio for individual cells depends on the isotope composition. Indicated in Fig. 14 are the cell filling fractions below . In this case, considering the nucleus in the shape of a droplet is still justified, see ref. . At , the nuclear matter distribution is described by pasta phases .
In the first part of this paper, we have performed nuclear structure calculations in Wigner-Seitz cells with uniform electron background. The calculations were carried out using the relativistic mean-field model with the NL3 parametrization. Physically motivated boundary conditions were introduced for the neutron single-particle states to describe unbound neutrons. We have studied the structural evolution of shapes along isotopic chains until very high neutron-rich systems and for fixed charge-to-baryon ratios. For superheavy isotopes, the generation of hollow shapes with a central deep in the proton distribution was found. We have also studied the single-particle states and mean-field potentials for different nuclei embedded in the neutron end electron background.
In the second part of the paper, we have replaced the constant electron background with realistic nonuniform distributions, found from the Poisson equation. The self-consistent calculations of nuclear and electric potentials were done within an iterative scheme, using the RMF description for the nucleons and the Thomas-Fermi approximation for the electrons. The differences between the electron chemical potentials for uniformly and nonuniformly distributed electrons were calculated for different atomic numbers and isotopes.
The whole nuclear chart was calculated for even numbers of protons between 8 and 240 and neutrons between to , both for realistic and uniform electron distributions. We have systematically compared the energies per baryon for these two cases and found the beta stability line as well as the neutron- and proton driplines. Additionally, we calculated the proton driplines for both cases and for the electron Fermi momentum of and , as well as for vacuum. At high electron densities, we have found a shift of several charge units in the proton dripline, caused by the realistic (non-uniform) electron density distributions. Despite the fact that all nuclear structure calculations in this paper were performed with RMF parametrization NL3, we believe that the main predictions will be similar for other popular versions of the RMF model.
This work is devoted to the memory of Prof. Walter Greiner, who was a great scientist, great teacher and great man. We started this work about 5 years ago after many fruitful discussions with him. We express our sincere gratitude for his continuous interest to this work and valuable advices. We also thank S. Schramm and U. Heinzmann for fruitful discussions. This work was supported in part by the Helmholz International Center for FAIR research.
-  C. J. Pethick and D. G. Ravenhall, “An Introduction to Matter at Subnuclear Densities,” in Neutron Stars: Theory and Observation, pp. 3–20, Dordrecht: Springer Netherlands, 1991.
-  N. K. Glendenning, Compact stars: nuclear physics, particle physics, and general relativity, vol. 242. 2000.
-  P. Haensel, A. Y. Potekhin, and D. G. Yakovlev, eds., Neutron Stars 1, vol. 326 of Astrophysics and Space Science Library. New York, NY: Springer New York, 2007.
-  G. Baym, C. Pethick, and P. Sutherland, “The Ground State of Matter at High Densities: Equation of State and Stellar Models,” The Astrophysical Journal, vol. 170, p. 299, dec 1971.
-  J. W. Negele and D. Vautherin, “Neutron star matter at sub-nuclear densities,” Nuclear Physics, Section A, vol. 207, no. 2, pp. 298–320, 1973.
-  P. Gögelein and H. Müther, “Nuclear matter in the crust of neutron stars,” Physical Review C, vol. 76, p. 024312, aug 2007.
-  F. Grill, J. Margueron, and N. Sandulescu, “Cluster structure of the inner crust of neutron stars in the Hartree-Fock-Bogoliubov approach,” Physical Review C, vol. 84, p. 065801, dec 2011.
-  H. S. Than, E. Khan, and N. Van Giai, “Wigner–Seitz cells in neutron star crust with finite-range interactions,” Journal of Physics G: Nuclear and Particle Physics, vol. 38, p. 025201, feb 2011.
-  B. K. Sharma, M. Centelles, X. Viñas, M. Baldo, and G. F. Burgio, “Unified equation of state for neutron stars on a microscopic basis,” Astronomy & Astrophysics, vol. 584, p. A103, dec 2015.
-  M. Baldo, E. E. Saperstein, and S. V. Tolokonnikov, “A realistic model of superfluidity in the neutron star inner crust,” The European Physical Journal A, vol. 32, pp. 97–108, apr 2007.
-  A. Pastore, S. Baroni, and C. Losa, “Superfluid properties of the inner crust of neutron stars,” Physical Review C, vol. 84, p. 065807, dec 2011.
-  S. Jin, A. Bulgac, K. Roche, and G. Wlazłowski, “Coordinate-space solver for superfluid many-fermion systems with the shifted conjugate-orthogonal conjugate-gradient method,” Physical Review C, vol. 95, p. 044302, apr 2017.
-  T. J. Bürvenich, I. N. Mishustin, and W. Greiner, “Nuclei embedded in an electron gas,” Physical Review C - Nuclear Physics, vol. 76, no. 3, 2007.
-  G. A. Lalazissis, J. König, and P. Ring, “New parametrization for the Lagrangian density of relativistic mean field theory,” Physical Review C, vol. 55, pp. 540–543, jan 1997.
-  C. Ebel, I. Mishustin, and W. Greiner, “Realistic electrostatic potentials in a neutron star crust,” Journal of Physics G: Nuclear and Particle Physics, vol. 42, p. 105201, oct 2015.
-  M. Bender, K. Rutz, P. G. Reinhard, and J. A. Maruhn, “Pairing gaps from nuclear mean-fieldmo dels,” The European Physical Journal A, vol. 8, pp. 59–75, jul 2000.
-  M. Baldo, E. Saperstein, and S. Tolokonnikov, “The role of the boundary conditions in the Wigner–Seitz approximation applied to the neutron star inner crust,” Nuclear Physics A, vol. 775, pp. 235–244, sep 2006.
-  N. Chamel, S. Naimi, E. Khan, and J. Margueron, “Validity of the Wigner-Seitz approximation in neutron star crust,” Physical Review C, vol. 75, p. 055806, may 2007.
-  S. Furusawa and I. Mishustin, “Self-consistent calculation of the nuclear composition in hot and dense stellar matter,” Physical Review C, vol. 95, p. 035802, mar 2017.
-  D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, “Structure of matter below nuclear saturation density,” Physical Review Letters, vol. 50, no. 26, pp. 2066–2069, 1983.
-  J. Lattimer, C. Pethick, D. Ravenhall, and D. Lamb, “Physical properties of hot, dense matter: The general case,” Nuclear Physics A, vol. 432, pp. 646–742, jan 1985.