Stability of spherical nuclei in the inner crust of neutron stars

Neutron stars are the densest objects in the Universe. In this paper we consider so-called inner crust - the layer, where neutron-excess nuclei are immersed into degenerate gas of electrons and sea of quasi-free neutrons. It was generally believed that spherical nuclei become unstable with respect to quadrupole deformations at high densities and here we consider this instability. Within perturbative approach we show that spherical nuclei with equilibrium number density are, in fact, stable with respect to infinitesimal quadrupole deformation. This is due to background of degenerate electrons and associated electrostatic potential which maintain stability of spherical nuclei. However, if the number of atomic nuclei per unit volume is much less than the equilibrium value, instability can arise. To avoid confusion we stress that our results are limited to infinitesimal deformations and do not guaranty strict thermodynamic stability of spherical nuclei. In particular, they does not exclude that substantially non-spherical nuclei (so-called pasta phase) represent thermodynamic equilibrium state of the densest layers of neutron star crust. Rather our results points that spherical nuclei can be metastable even if they are not energetically favourable and the timescale of transformation of spherical nuclei to the pasta phases should be estimated subsequently.


Introduction
The inner crust of neutron stars extends from the density ρ drip ≈ 4.3 × 10 11 g/cm 3 to ∼ 10 14 g/cm 3 [1,2]. It consists of fully ionised atomic nuclei immersed into background of quasi-free neutrons and relativistic degenerate electron gas. Atomic nuclei have large neutron excess because of the high chemical potential of the electrons. They have a spherical shape in most of the inner crust, but in the deepest layers energetically favourable nuclei configuration can become substantially non-spherical (cylinders, planes and inverse configurations; so-called pasta phases) [1][2][3][4][5][6] (see also [7][8][9][10] for the most recent progress). 1 Following [13], the density region for spherical nuclei was generally assumed to be limited by instability of spherical nuclei with respect to quadrupole deformations (see, e.g., [1,2,14]). Namely, applying the instability criterion derived by Bohr and Wheeler [15] for isolated nucleus, the spherical nuclei were predicted to become absolutely unstable (even for infinitesimal deformation) when the ratio of the nucleus volume to the Wigner-Seitz cell volume (filling factor) u reaches a value of 1/8 = 0.125 [13]. Indeed, recent extended Thomas-Fermi calculations of Ref. [8] reports a transition from spherical to cylindrical nuclei at filling factor close to 0.125 (see their table XII and respective discussion). However, calculations within the compressible liquid drop model (CLDM) typically predict spherical nuclei to be energetically favourable up to the larger filling factor ∼ 0.2 (see, e.g., [4,11] and Fig. 1 for numerical illustration). Thus, at least within CLDM approach, there is a contradiction between numerical results and predicted instability: spherical nuclei are predicted to be thermodynamically stable at the region there they supposed to be absolutely unstable. Obviously, it have two (not mutually exclusive) 1 Some works (e.g., [11,12]) predict pasta phases to be absent in neutron stars. solutions: A) the true thermodynamic equilibrium for filling factors 0.125-0.2 correspond to complex nuclear structures (e.g., [9]), which stays beyond the scope of most works based on CLDM approach; B) the instability of spherical nuclei is suppressed in the inner crust. The option B), enhanced stability of spherical nuclei, were suggested by a number of authors (e.g., [11,20]). The strongest point in support of this solution was suggested in Ref. [21], which, in particular, investigated the stability of spherical nuclei with respect to infinitesimal deformations within the Wigner-Seitz cell approximation. 2 Authors concluded that spherical nuclei are stable at arbitrary filling factor, if the number density of nuclei corresponds to the optimal value. 3 However, in our opinion, results of Ref. [21] are based on inaccurate boundary conditions. Namely, to ensure that the electric field flux over the cell boundary is zero, authors of Ref. [21] impose the Nuemann boundary condition, i.e. demand that the normal component of the electric field is zero at each point of the cell boundary. The latter seems unreasonable: thanks to the zero charge of the cell and Gauss theorem, the electric field flux over the cell boundary vanishes automatically if electron and nucleus (proton) contributions are both correctly incorporated. Imposing of specific boundary conditions, in particular, the Neumann boundary condition, is equal to assumption that the charges outside the cell rearrange themselves to a distribution which is required to ensure imposed boundary condition. Note, that the required distribution of outside charges depends on deformation of nucleus, making such rearranging unnatural from our point of view. Here we correct this problem. In addition, our consideration naturally takes into account the neutron skin, which is essential ingredient of two phase system boundary thermodynamics (e.g., [24,25]), however, account of this effect do not change the results.
To avoid confusion, let us stress that we analyse only stability with respect to infinitesimal deformations. Clearly, it does not sufficient to guaranty absolute thermodynamic stability. Indeed, as shown by Hashimoto et al. [4], the spherical nuclei can not be energetically favourable at too large filling factors and thus spherical nuclei can be treated only as metastable in strict thermodynamic sense.

Calculations
In order to explain the nature of the above mentioned contradictions we checked validity of the option B): suppression of spherical nuclei instability in inner crust. We applied CLDM of Ref. [26], in which the nucleus is surrounded by quasi-free neutrons, being located at the center of a spherical Wigner-Seitz cell, which is uniformly filled with electrons (according to the quasi-neutrality requirement the total charge of the cell is zero). This model naturally incorporates neutron skin effects (see supplementary material to Ref. [26] for details).
Following Bohr and Wheeler [15], we considered the change of energy when the nucleus is deformed from a sphere to an spheroid with semi-axes R(1 + ε) and R/ √ 1 + ε. Here R is the radius of the spherical nucleus, ε is infinitesimal deformation parameter. The Wigner-Seitz cell assumed to stay spherical. Within CLDM of Ref. [26] the change of the cell energy can be calculated analytically up to the second order in ε (see Appendix for the details): Here σ is the surface tension, R is the radius of the nucleus before deformation, Ze is the nucleus charge, e is the elementary charge. Eq. (1) is similar to Eq. (9) of Ref. [15], but does not coincide exactly (even for ε 2 terms, considered here) due to presence of electron background, which induces electrostatic potential and modifies the Coulomb energy change associated with nuclei deformation. The difference is ∝ u. It agrees with results of Ref. [22], as they were cited in [13]. However, as we demonstrate below, this term can not be neglected because it leads to suppression of the instability. It is worth to stress that Eq. (1) holds true in exactly the same form for a simplified CLDM, which neglects neutron skin effects (see Appendix for derivation details), so our results are equally applicable for this widely applied type of CLDMs.

Equilibrium inner crust
The virial theorem [23] (see also supplementary material in [26] for derivation with accurate account of neutron adsorption effects) provide coupling for Coulomb energy and surface terms, if crust composition corresponds to the equilibrium (the number of atomic nuclei per unit volume is optimal) : Substituting 4πσR 2 , given by (2), into Eq. (1) we got the energy change associated with infinitesimal nuclei deformation in equilibrium crust: It is easy to check that at any value of the filling factor u ∈ (0, 1) the energy change δE remains positive. Thus, the nuclei in the crust with equilibrium composition remain stable with respect to infinitesimal quadrupole deformations at any values of the filling factor. It is worth to stress that this result is derived analytically and stays the same for CLDMs which includes and neglects neutron skin effects. In particular, it does not depend on the choice of nuclear physical model required to specify numerical parameters of CLDM (e.g. the bulk energy and surface tension). Note that the equilibrium composition is believed to be a good model for nonaccreting neutron stars [1] (see however [27,28]).

Non-equilibrium inner crust
In the more general case the number of atomic nuclei per unit volume can differ from equilibrium value. In this case, instead of (2) we can write more general expression [26]: Here µ N is the the chemical potential of the nucleus, which describes the change in energy when one nucleus is created from nucleons, which are already present in the substance (for equilibrium crust µ N = 0). Combining equations (1) and (4), we obtain the energy change associated with deformation of nucleus in non-equilibrium crust For accreting neutron stars, there is an excess of nuclei in the crust [26]. It leads to µ N > 0 and makes the spherical nuclei even more stable with respect to quadrupole deformations. 4 However, if µ N < 0, i.e., when the number density of nuclei in the stellar matter is less than the equilibrium value, the nuclei may become unstable with respect to quadrupole deformations and this instability likely leads to fission.
For numerical illustration of possible instability we perform calculations including neutron skin (adsorption) effects within CLDM of Ref. [26] and apply SLy4 nucleon-nucleon potential [16,17]. We parametrize matter by baryon number density n b and parameter C = µ N /µ n , which stays constant in the inner crust thanks to diffusion/hydrostatic equilibrium of quasifree neutrons [26]. The results are shown in Fig. 2 for several values of C.  [26]. SLy4 potential [16,17] is applied.

Discussion, results and conclusions
Within CLDM we demonstrate that spherical nuclei are stable with respect to infinitesimal quadrupole deformations, if their number density correspond to the equilibrium value. The suppression of the instability, in comparison with isolated nuclei, is due to the fact that nuclei in the inner crust are immersed into background of degenerate electrons. The electron charge density is comparable with the charge density of the nucleus and induces electrostatic potential, which supports spherical shape of nuclei.
For non-equilibrium crust, when the number of nuclei per unit volume is less than the equilibrium number, the instability with respect to infinitesimal quadrupole deformations can appear. This results seems to be natural: fission, likely caused by instability, leads to an increase in the nuclei number, driving composition closer to the equilibrium. According to our calculations for SLy4 potential [16,17], nuclei number density should be lower than the equilibrium value for a factor of 2.2-2.4 to ensure the instability at n b > 0.047 fm −3 . If nuclei number density is larger than in the equilibrium at the same baryon density (e.g., in the accreted crust [26]), the spherical nuclei are stable with respect to infinitesimal deformation. However, it does not exclude other types of instability which are not associated with deformation of nuclei, for example, the instability considered in [26].
Qualitatively similar results were obtained in [21], but as we point in introduction, this work appeal to artificial boundary conditions. It leads to quantitative differences with our results.
It is worth to stress that our analysis is perturbative, and thus limited to infinitesimal deformations. Obviously, it can not guaranty absolute thermodynamic stability of spherical nuclei. And indeed, the cylindrical and other pasta phases shown to become more energetically favourable at large filling factors, u 0.2 [4]. Combination of stability of spherical nuclei for infinitesimal deformations with instability in strict thermodynamic sense for u 0.2 suggest that in this conditions spherical nuclei should be, in fact, metastable (i.e. correspond to local energy minimum, which differs from the global minimum). It opens an interesting task to estimate the transition timescale from spherical to nonspherical shapes, however, we leave any estimates for this timescale beyond of this work. However, we point that u = 1/8 criterion should not be applied to put an upper bound for the density region of spherical nuclei in equilibrium crust because it was suggested on the base of simplified consideration of fission instability, which is not supported by more detailed analysis, which was performed in our work.
Concluding, we should warn the reader than all our results are obtained within spherical Wigner-Seitz cell approximation, which is likely inaccurate and should be elaborated for very large filling factor u ∼ 1. However, in realistic models of crust the filling factor for spherical nuclei is u 0.2 and spherical Wigner-Seitz cell seems a reasonable approximation, but we don't check the latter statement straightforwardly. We also warn that we neglect the curvature corrections to the surface tension, which shown to be important for thermodynamically determined boundaries of the pasta layers (e.g., [29]), however, we don't expect that it can affect our results qualitatively.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: At the appendix we derive Eq. (1) which describes the change of the cell energy associated with deformation of nuclei. The derivations are provided for two versions of CLDM, one which neglects neutron skin effects and model of Ref. [26], which incorporates these effects thermodynamically consistently. The curvature corrections are neglected in both cases. In the first subsection A.1 we derive the Coulomb energy of a cell with deformed nuclei, which is required for derivation of Eq. (1).

Appendix A.1. Coulomb energy of the cell with deformed nucleus
We start from consideration of the Coulomb energy of a cell with deformed nucleus in the center. It is essentially electrostatic problem: calculation of the Coulomb energy of uniformly negatively charged sphere with radius R c = R/u 1/3 and uniformly positively charged spheroid (with semi-axes R(1 + ε) and R/ √ 1 + ε), located at the center of the sphere. The net charge of this system is zero. We neglect terms ∝ ε 3 .
We begin with the Poisson's equation for Wigner-Seitz cell where where r and θ are radial distance and the polar angle of spherical coordinate system (r = 0 is the center of the cell), Θ(x) is the Heaviside step function. ρ pi = Ze/(4πR 3 /3) is the proton charge density inside nucleus, ρ e = −uρ pi is the electron charge density (the cell is electrically neutral). The protons are located within a spheroid, with boundary given by It is worth to stress, that ρ p does not depend on ε because we consider deformation which do not affect nucleus volume. The solution to Eq. (A1) can be presented as a sum of the proton potential ϕ p and the electron potential ϕ e , while the total Coulomb energy can be presented as , and E e−e C are proton-proton, electron-proton, and electronelectron contributions, respectively: Here we take into account that proton density is zero outside the spheroid. The explicit form of electron potential inside cell is well known while the proton potential inside the nucleus (r < R sp (θ)) can be written in form: The omitted terms contributes to the Coulomb energy only at order ε 3 or higher, which is not considered here. The proton potential outside nucleus is not required to calculate Coulomb energy (see Eqs. A3-A5) and thus not shown here. The integrals in Eqs. (A3)-(A5) can be calculated analytically. Up to the second order in ε they are: For ε = 0 they agree with well known expression for the Coulomb energy of the cell with spherical nucleus (e.g., [23]).

Appendix A.2. Calculation of the Energy change neglecting neutron skin
Within CLDM, which neglects neutron skin the surface tension is typically assumed to be function of x i -the ratio of proton number density to the total baryon number density inside the nucleus. To derive Eq. (1) we consider difference of the cell energies between two configurations: (a) WS cell with spherical nucleus and (b) WS cell with deformed nucleus.
We assume that the proton and neutron number densities inside nuclei as well as neutron number density outside nucleus are the same in both configurations. The volume of nucleus is also unmodified by deformation. In this case the energy change associated with deformation, δE = E b − E a , contain only surface and Coulomb terms, while the bulk terms are cancel out.
The surface tension is the same in both configurations (x i is not modified) and the surface energy change is given by difference between the surface areas of the spheroid and the sphere. Thus, change of the surface energy is The change of the Coulomb energy can be calculated using results of section A.1, which leads to The net energy change δE = δE C + δE s , thus summing up Eq. (A12) and (A11), we obtain Eq. (1).

Appendix A.3. Calculation of the energy change including neutron skin
Derivation of Eq. (1) for CLDM of Ref. [26], which accounts for neutron skin (adsorption) effects, is more complicated. It is due to the fact that nuclei deformation changes of nuclei surface area and thus amount of neutrons adsorbed on it. Thus, one can not assume that neutron number densities and nucleus volume are not modified by deformation because it will lead to variation of total number of neutrons (and thus net baryon density).
To derive Eq. (1) we consider an extended version CLDM of Ref. [26], which is allowed for deformed nuclei. Namely, nucleus deformation parameter ε is added as additional independent CLDM variable. It allows to write down the expression for the energy density , which differs from Eq. (2) of the supplementary material to Ref. [26] only by the Coulomb and surface energy terms: Here bulk (n n , n p ) is the energy density of bulk nuclear matter at respective neutron and proton number density (n n and n p ), n ni , n pi are the neutron and proton number density inside nucleus, n no is the neutron number density outside nucleus (we assume that proton drip does not take place), and e e (n e ) is the energy density of degenerate electrons at electron number density n e = u n pi . The cell volume V c = 4πR 3 /(3u). According to the results of section A.1, the Coulomb energy is where f (u) = 1 − 3 u 1/3 /2 + u/2. The surface energy where µ ns and ν s are chemical potential and surface density of the adsorbed neutrons. The nuclei surface area A is given by the surface area of the spheroid The thermodynamic consistency requires dσ dν s = −ν s dµ ns dν s . (A17) Thus σ and µ ns can be treated as functions of ν s . Following [26] let us introduce auxiliary variables, which simplifies subsequent analysis Here N s = Aν s is the number of neutrons, adsorbed to the nucleus.
Here we analyse stability with respect to infinitesimal deformation of nuclei, described by parameter ε. To do so we check for two conditions: (a) ε = 0 is an extremum of the energy density and (b) the extremum at ε = 0 is a (local) minimum. While doing this analysis we assume that net baryon number density ns , and u, generally, can vary with variation of ε. As the first step, we need to specify equlibrium values of internal CLDM parameters, which are given by the condition that the partial derivatives with respect to each of internal parameter is zero. As long as we are interested on these parameters at ε = 0, i.e. for spherical nuclei, respective equations are exactly the same as in [26]. In particular, if n N differs from the equilibrium value, Eq. (4) holds true.
Let us note that the energy density explicitly depends only on ε 2 , thus the partial derivative ∂ /∂ε ∝ ε, thus ε = 0 is indeed an extremum.
To check is the extremum at ε = 0 minimum or maximum we, firstly, note that the mixed partial derivatives ∂ 2 /∂ε∂p ∝ ε, and thus they are zero, if calculated at ε = 0 (here p is arbitrary parameter of the CLDM model, except ε). As a result, the variation of energy density associated with infinitesimal ε is The bulk terms in the energy density do not depend on ε explicitly and only Coulomb and surface energy contributes to ∂ 2 /∂ε 2 . These derivatives can be easily calculated, 5 leading to δ = n N 8π 5 R 2 σ + 3 5 Z 2 e 2 R u 2 − 1 5 ε 2 . (A25) 5 As long as depends on ε only via ε 2 , ∂ 2 /∂ε 2 ε=0 = 2 ∂ /∂ ε 2 ε 2 =0 . While calculating derivative of the surface energy, one should take in mind that it is calculated at fixed n (tot) ns and thus fixed total amount of adsorbed neutrons. In this case δE s = σδA due to Eq. A17.