text:
# Structure of the inner crust of neutron stars: crystal lattice or disordered phase ?

## I INTRODUCTION

The outer layer of a neutron star has a complex structure which depends strongly on the nuclear density. In the inner crust of the star, due to the high density and pressure, a large fraction of neutrons occupies unbound states. Nuclei which are still present are therefore immersed in a neutron gas. The structure of this part of the star has been the subject of a considerable theoretical effort . The main conclusion emerging from these studies is that the nuclei form exotic phases, stabilized by the Coulomb interaction, and characterized by various crystal lattice structures.

The theoretical predictions on the structure and stability of the nuclear phases are based on the minimization of the energy of the neutron-proton-electron (npe) matter within a Wigner-Seitz cell, whose radius is treated as a variational parameter. An agreement has been reached on the existence of five phases formed in the region where the nucleon density varies from $`0.03`$ to $`0.1`$ fm-3. The following chain of phase transitions is predicted as the density increases: spherical nuclei $``$ rods $``$ slabs $``$ tubes $``$ bubbles $``$ uniform matter. The physical mechanism responsible for the phase transition pattern results from the interplay between the Coulomb and surface energies. However, the differences between the energy density of the various phases are only of the order of a few $`keV/fm^3`$, for a wide range of densities . The precise transition pattern is therefore dependent on the ingredients of the theoretical models. For example, recent calculations based on the self-consistent Hartree-Fock (HF) method have predicted that the phase composed of spherical nuclei is stable up to a density around $`0.077`$ fm-3, above which the uniform matter is energetically favored . Other phases appear at higher energies and form excited configurations.

In the present paper we focus on an effect which has been overlooked in previous studies but should strongly influence the phase transition pattern. In all to date approaches, shell effects have been either completely neglected, as in models based on a liquid drop formula or in Thomas-Fermi models, or have been taken into account only for bound nucleons. On the other hand, it has been recently shown in the semiclassical approximation, based on the Gutzwiller trace formula, that the shell effects in inhomogeneous nuclear matter are quite strong and yield corrections to the energy density of the order of several $`keV/fm^3`$. This value is comparable to the differences between the liquid drop energy densities predicted for the different nuclear phases . Moreover, the shell energy oscillates as a function of density and depends strongly on the spatial order in the system. In particular, its amplitude depends both on the nuclear geometry and on the lattice type.

The shell energy of npe matter in the neutron star crust has two origins. The first one is associated with bound nucleons and has been shown to be small by Oyamatsu and Yamada in Ref. . The second one is related to unbound neutrons which may scatter on inhomogeneities and form resonant states. The situation is somewhat similar as for the Casimir effect in quantum field theory and condensed matter (see e.g. Refs.  and references therein), where, in some systems, a fluctuation induced interaction leads to a correction to the total energy. In the same way, the scattering of unbound neutrons lead to an effective interaction between nuclei immersed in the neutron environment. This effect can be termed the fermionic Casimir effect and has been studied for simple geometries in Refs. . For impurities that are small compared to the neutron Fermi wavelength, the generated interaction is mainly due to s-wave scattering. It gives then rise to the Ruderman-Kittel correction which is small compared to the Fermi energy of neutrons .

However, the s-wave scattering limit is not valid in the neutron star crust, since the neutron Fermi momentum $`k_F^n`$ is of the order of $`1fm^1`$ and thus $`k_F^nR>1`$, where $`R`$ is the size of an inhomogeneity. Consequently the contribution of higher partial waves to the scattering process cannot be neglected, and gives rise to much larger energy correction, as for large objects more of the incident wave will be reflected . The leading shell energy contribution, associated with the neutron scattering, to the total energy can be determined in the semiclassical approximation for various geometries of nuclear shapes (see appendix B). It has been shown to yield energy corrections that influence the structure of inhomogeneous nuclear matter . It is also obvious that since the gross shell structure is determined by the shortest periodic orbits in the system, the shell energy is a sensitive function of the crystal lattice type, as well as of the shape of inhomogeneities . Unfortunately the semiclassical approach has drawbacks. First, it assumes that the obstacles are impenetrable scatterers which may overestimate the amplitude of the shell energy. Second, the method lacks of the mutual interplay between the shell energy and the liquid drop energy. Such a coupling is quite important since it allows a part of the shell energy to be “absorbed” into deformation. Hence a microscopic treatment of the problem is needed where both effects, i.e. the one coming from the liquid drop energy and the shell energy term are treated on the same footing.

## II THE METHOD

The ground and excited (isomeric) states of npe matter have been determined by the Skyrme-Hartree-Fock method. We have solved the HF equation by discretization in coordinate space, within a rectangular box and with periodic boundary conditions imposed on the nucleon wave functions . This technique has the advantage that it allows to describe any kind of shapes of the nuclear density distribution with the same accuracy, and to introduce naturally periodic boundary conditions. The code that was previously developed to study nuclei with an ellipsoidal symmetry has been modified to accommodate the periodic solutions.

The high efficiency of the method can be attributed to the iterative procedure used to solve the HF equations, the imaginary time step method . One of its features is that only those single-particle wave functions which contribute to the mean field are calculated. It is particularly important in our case, since the densities which we are dealing with here require to consider up to $`1000`$ nucleons per box. In order to take them into account properly, we have used $`600`$ wave functions. Due to Kramers degeneracy, they correspond to $`1200`$ single-particle states and thus form a sufficiently large Hilbert space for our purposes.

As a starting point of the HF iterations, Fermi gas wave functions have been used, which turned out to provide a good first guess of the true ground state. The criterion for the termination of the iteration process is the total energy difference obtained after a fixed number of imaginary time steps. In our calculations we terminated the iteration procedure when the differences between the total energy calculated at two successive iterations did not exceed $`1keV`$. For stable configurations this value was of an order of magnitude lower.

An important parameter of a mesh is the spacing between the mesh points. It must be small enough to lead to accurate results and large enough to accommodate the very large number of wave functions needed to describe the neutron gas. We have chosen as a compromise between these two requirements a value of the mesh size equal to $`\mathrm{\Delta }x=1.3fm`$.

This choice of $`\mathrm{\Delta }x`$ fixes a bound to the maximum momentum of a nucleon that can be represented by the discretized wave function. The highest momentum is indeed equal to $`k_{crit}=\pi /\mathrm{\Delta }x2.42fm^1`$, and corresponds to the maximum kinetic energy $`ϵ_{crit}={\displaystyle \frac{\mathrm{}^2\pi ^2}{2m\mathrm{\Delta }x^2}}121MeV`$.

There are two main sources of errors caused by the discretization. One of them is related to unbound neutrons which may occupy high energy states and thus can not be properly described by a discretized wave function. However even at the highest neutron density that we are interested in here, i.e. $`\rho _n0.1`$ fm-3, the neutron Fermi momentum is equal to $`k_F^n1.44`$ fm-1 which is far below the critical value. A second source of errors is associated with nucleons that are bound in a nucleus. Typically, a nuclear potential at densities $`\rho >0.03`$ fm-3 has a depth which does not exceed $`20`$ MeV and a radius lower than $`10`$ fm . In a Fermi gas approximation for bound nucleons, the momentum of the last bound nucleon is around $`1`$ fm-1, well below the maximal value that can be represented on the mesh. Let us note that the Fermi gas approximation gives an upper bound of the maximum momentum since the wave function of a bound state has a long tail due to the large diffuseness of a nuclear potential immersed in a neutron gas.

To decrease the numerical task, we have assumed that the nuclear density distribution is symmetrical with respect to three planes, allowing however the existence of triaxial shapes. Parity breaking configurations are thus excluded.

Different phases may coexist for a given nuclear density. To determine them and to study their stability requires the construction of adiabatic paths between phases. This has been done by introducing a constraint on the quadrupole moment of the proton density distribution in the Skyrme HF equations. In order to compare energies of configurations differing by the proton quadrupole moment, a careful treatment of the Coulomb interaction, both between protons in different nuclei of the lattice, as well as between protons and electrons, is needed. It has lead us to solve the Poisson equation for the electric potential within a box with periodic boundary conditions. This approach is clearly free from limitations of the usual Wigner-Seitz approximation, where the influence of the lattice type on the Coulomb energy is not taken into account. The electrons are assumed to form a uniform relativistic gas where the screening effects are neglected. The total charge of the box is zero, which gives a relation between the electron and proton densities.

Summarizing, we solve the Hartree-Fock equations for the npe matter with contributions in the energy functional coming from the nuclei, the neutron gas and the electrons. In such a calculation, the liquid drop and the shell energy parts of the energy are automatically and self-consistently included. The total energy of the system can be expressed in the form:

$$E_{tot}(\rho _p,\rho _n,\rho _e)=E_{Skyrme}(\rho _p,\rho _n)+E_{Coul}(\rho _p,\rho _e)+E_{el}(\rho _e),$$
(1)
where $`E_{Skyrme}`$ comes from the Skyrme energy functional , $`E_{Coul}`$ is the Coulomb energy term, and $`E_{el}`$ denotes the electron kinetic energy term. The densities $`\rho _n,\rho _p,\rho _e`$ denote the average densities of neutrons, protons and electrons, respectively. They are related to the density distribution through the equation:

$$\rho _i=\frac{1}{V}_V\rho _i(𝐫)d^3r,\text{ i=n,p,e},$$
(2)
where $`V`$ is the volume of the cell. The charge neutrality condition implies that $`\rho _e=\rho _p`$. The total nuclear density is given by $`\rho =\rho _n+\rho _p`$. In the calculations we have used the SLy4 Skyrme force. This recent parametrization  gives a good description of neutron matter at low densities . For a few test cases, calculations have also been performed with the SLy7 force which does not lead to any qualitative differences.

The shell energy which is implicitly included in the expression (1) can be decomposed into two parts. The first one, denoted by $`E_{shell}^{in}`$, is determined by the density of neutrons $`\rho _n^{in}`$ and protons $`\rho _p^{in}`$ inside the nucleus i.e.:

$$\rho _i^{in}=\frac{1}{V_{nucl}}_{V_{nucl}}\rho _i(𝐫)d^3r,\text{ i=n,p},$$
(3)
where $`V_{nucl}`$ is the volume of the nucleus. The second, $`E_{shell}^{out}`$, is a function of unbound neutrons $`\rho _n^{out}`$ and protons $`\rho _p^{out}`$, where

$$\rho _i^{out}=\frac{1}{VV_{nucl}}_{VV_{nucl}}\rho _i(𝐫)d^3r,\text{ i=n,p}.$$
(4)
Since $`\rho _p^{out}`$ is practically zero for $`\rho <0.08`$ fm-3 , the contribution to the energy due to unbound protons is marginal.

It should be emphasized that the decomposition of the density into an inner and an outer part is somewhat arbitrary. There are various methods of defining these quantities that have been used in liquid drop based approaches . In the HF method this decomposition is not necessary to calculate the total energy of the system. Nevertheless, in order to understand its behavior as a function of physical parameters such as nuclear density, proton fraction, lattice constant, etc., it is important to identify and analyze the origin of the various contributions to the total energy.

## III RESULTS

The HF equations are solved for a fixed total density $`\rho `$, box size $`d`$ and with a constraint imposed on the proton quadrupole moment. In this way, it is possible to determine not only the ground state but also excited configurations, differing from the ground state by, for example, the lattice type. Since we solve the problem in a cubic box with three symmetry planes, several lattice geometries can be generated like e.g. simple cubic crystal (scc), face centered crystal (fcc), or body centered crystal (bcc) .

In order to verify that the absolute minimum of the functional (1) has been found, we have performed a set of constrained calculations corresponding to various values of the proton quadrupole moment for each value of the total nuclear density. It should be emphasized that we did not perform the total energy minimization with respect to the parameters of the box which were kept fixed at values corresponding to those obtained in the reference .

In the figure 1 we show the ground state energy density as a function of the proton fraction: $`Z/A=\rho _p/\rho `$ ($`Z=V\rho _p`$ and $`A=V\rho `$, where $`V=d^3`$ is the volume of the box) for four values of the nucleon density. The calculations indicate that for the density range we consider here, the optimal $`Z/A`$ ratio varies between $`0.03`$ and $`0.04`$.

For a total density $`\rho =0.055`$ fm-3, the energy density is plotted in Fig. 2 (middle part of the left panel) as a function of the proton quadrupole moment $`Q_p=Q_{20}^p`$. Most studies predict that, at this density, the nuclear spherical phase is energetically favored (see e.g.  and references therein). One can see that in our self-consistent calculations, the energy minimum appears for a non-zero quadrupole moment. In order to identify the shape of the nuclei, as well as the lattice type, we have plotted the neutron and proton density distributions integrated along the $`z`$ and $`y`$ axes:

$`\rho _i^z(x,y)`$ $`=`$ $`{\displaystyle _0^d}\rho _i(x,y,z)𝑑z,`$ (5)
$`\rho _i^y(x,z)`$ $`=`$ $`{\displaystyle _0^d}\rho _i(x,y,z)𝑑y,`$ (6)
where $`i=n,p`$. In the lower and upper parts of the left panel of Fig. 2, the integrated density distributions are plotted for both protons and neutrons. We have shown the density distributions for the configurations indicated by arrows, which correspond to the spherical and the most elongated nuclear shapes. The density distributions in the energy minimum (not plotted) have a pattern very similar to the spherical one. They correspond to a scc lattice of deformed nuclei immersed in a neutron gas.

An excited minimum appears for a proton quadrupole moment of about $`2500fm^2`$. It is associated with a rod-like (“spaghetti”) phase (upper part of the left panel of Fig. 2). The barrier between these minima is of the order of $`0.5keV/fm^3`$. One can see that contrary to previous calculations, the spherical phase is no longer the most favored one and appears at an excitation energy of about $`1keV/fm^3`$, even above the “spaghetti” phase. When keeping the proton number fixed and adding neutrons to the system, the energy density pattern slightly changes. The right panel of Fig. 2 corresponds to the situation when the number of neutrons have been increased to reach the total density $`\rho =0.057`$ fm-3. As shown in the right panel of Fig. 2, the position of the minimum is not altered. It suggests that the deformation of the ground state is mainly related to the bound nucleons. On the other hand, there is a substantial change of the position of the spherical phase which drops down in energy and is now below the “spaghetti” phase.

The Fig. 2 shows that the change in the neutron density does not influence qualitatively the structure of the ground state since the energy minimum stays approximately at the same proton quadrupole deformation. However it is not always the case as one can see from Fig. 3, where the energy density has been plotted as a function of proton quadrupole moment for larger total densities. The left panel shows the results for $`\rho =0.075`$ fm-3, and the right one for $`\rho =0.079`$ fm-3. The two almost degenerate minima visible in the left panel differ by the lattice type only, and both correspond to nuclei immersed in a neutron gas. The spherical configuration corresponds to a scc lattice and the deformed one to a bcc lattice. There is also an excited configuration associated with a slab-like (“lasagna”) phase. Still increasing the number of neutrons, but at fixed proton number, one sees on the right panel of Fig. 3, that the relative position of the minima is not stable with respect to the neutron density. The spherical phase is no longer energetically favored and the bcc configuration merges with the slab phase, forming the lowest energy configuration in the form of “rippled” nuclear slabs.

The question arises whether these changes are related to the shell effects or to the liquid drop part of the energy. In order to better visualize this effect, we have plotted (see Fig.4) the energy differences between different configurations:

* the spherical nuclear phase ($`Q_p=0`$ fm2),
* the “spaghetti” phase ($`Q_p=2400`$ fm2),
* the “lasagna” phase ($`Q_p=2600`$ fm2).

The solid curve represents the quantity $`\mathrm{\Delta }E/V=(E_{spherical}E_{spaghetti})/V`$ while the dotted curve corresponds to $`\mathrm{\Delta }E/V=(E_{spherical}E_{lasagna})/V`$. One notices that these differences oscillate rapidly as a function of the total density. It seems therefore unlikely that the liquid drop part of the energy is responsible for such rapid fluctuations. On the other hand, the leading shell energy contribution to the total energy, for two obstacles located at the distance $`d`$ and immersed in the neutron gas, can be determined in the semiclassical approximation (see appendix B):

$$E_{shell}^{out}\frac{\mathrm{}^2L^iR^{2i}}{8m_n}\left(\frac{3}{\pi }\right)^{\frac{2+i}{6}}\frac{(\rho _n^{out})^{\frac{2+i}{6}}}{d^{\frac{6i}{2}}}\mathrm{cos}\left(2k_F^ndi\frac{\pi }{4}\right),$$
(8)
where $`m_n`$ is the neutron mass and $`i=0,1,2`$ stands for two spherical, cylindrical and planar obstacles, respectively (we assume that rods and slabs are parallel to each other). In the above equation, $`L`$ defines the length of the obstacle, and $`R`$ is its radius (in the case of a slab it is defined as half of its width). This expression exhibits oscillations, where the gross staggering pattern comes from the shortest periodic orbit of length $`d`$ between nearest nuclei. This can be used to estimate qualitatively how the difference between the shell energies of the two minima varies as a function of the neutron density.

Assuming that the distance between the nearest nuclei is given by $`r`$, the shell energy contribution due to unbound neutrons has the form:

$$E_{shell}^{out}\mathrm{cos}(2k_F^nr).$$
(9)
The quantity plotted in Fig. 4 represents the difference between the energy density of two phases, and if it is associated with shell effects, its behavior should be expressible as a sum of cosine (or sine) functions. In particular, the leading contribution with the longest period comes from the single repetition of the shortest periodic orbits:

$$\mathrm{\Delta }E_{shell}^{out}A_i\mathrm{cos}(2k_F^{ni}r_i)+A_j\mathrm{cos}(2k_F^{nj}r_j)\text{ }ij,$$
(10)
where $`i,j=0,1,2`$ for spherical nuclear phase, the “spaghetti” phase, and the “lasagna” phase, respectively. The coefficient $`A_i`$ is related to the stability of the orbit, $`r_i`$ is the distance between nuclear inhomogeneities, and $`k_F^{ni}`$ stands for the Fermi momentum in the phase $`i`$. One can infer from this expression that the interference between the contributions of different orbits is at the origin of the fluctuations of $`\mathrm{\Delta }E_{shell}^{out}`$ as a function of the density.

Indeed, taking $`k_F^{ni}=k_F^n(3\pi ^2\rho _n^{out})^{1/3}`$ and assuming that $`A_iA_j=A`$ we get:

$$\mathrm{\Delta }E_{shell}^{out}2A\mathrm{sin}(2k_F^nr_{ij})\mathrm{sin}(k_F^n\mathrm{\Delta }r_{ij}),$$
(11)
where $`r_{ij}=\frac{1}{2}(r_i+r_j)`$ and $`\mathrm{\Delta }r_{ij}=(r_ir_j)`$. Hence the shell energy difference oscillates as a function of the Fermi momentum with a period $`\mathrm{\Delta }k_F^n={\displaystyle \frac{\pi }{r_{ij}}}`$, or equivalently, as a function of the neutron Fermi energy, with the period $`\mathrm{\Delta }ϵ_F^n={\displaystyle \frac{\mathrm{}^2\pi ^2}{2m_nr_{ij}^2}}`$.

This quantity can be estimated from our HF calculations. Let us define the radius of the nuclear inhomogeneity as the distance between the center of a nucleus and the point where the neutron density decreases to the value $`\frac{1}{2}(\rho _n^{out}+\rho _n^{in})`$. This gives a radius $`R_06.23fm`$ for the spherical and $`R_13.86fm`$ for the “spaghetti” phases, at a density $`\rho 0.055`$ $`fm^3`$. Since the length of the lattice $`d`$ is in this case equal to $`26fm`$, the distance between the two nearest spherical-like and spaghetti-like objects can be estimated as $`r_0=d2R_013.54fm`$, and $`r_1=d2R_118.28fm`$, respectively. This gives a period of oscillation equal to $`\mathrm{\Delta }ϵ_F^n0.8`$ MeV. On the other hand, the HF calculations predict the neutron Fermi energies to be equal to:

$`ϵ_F^n(Q_p=0,\rho =0.063)`$ $`=`$ $`12.043MeV,`$ (12)
$`ϵ_F^n(Q_p=0,\rho =0.053)`$ $`=`$ $`10.74MeV,`$ (13)
$`ϵ_F^n(Q_p=2400,\rho =0.063)`$ $`=`$ $`11.931MeV,`$ (14)
$`ϵ_F^n(Q_p=2400,\rho =0.053)`$ $`=`$ $`10.607MeV,`$ (15)
which give the periods:

$`\mathrm{\Delta }ϵ_F^n(Q_p=0)`$ $`=`$ $`1.303MeV,`$ (16)
$`\mathrm{\Delta }ϵ_F^n(Q_p=2400)`$ $`=`$ $`1.324MeV.`$ (17)
They have to be compared with the $`0.8`$ MeV obtained on the basis of the semiclassical analysis. The difference between these values is mainly related to the fact that in the semiclassical approximation we have neglected both the coupling of the shell effects to the liquid drop terms, and higher order corrections related to other orbits (including the multiple repetitions of the shortest orbits). Nevertheless, taking into account how drastic these approximations are, the agreement is surprisingly good. It is even better for the case of the energy density difference between the spherical and the slab phases (dotted curve in the Fig. 4). The radii are then: $`R_05.9fm`$, and $`R_23.59fm`$ at a density $`\rho 0.075`$ $`fm^3`$, with the “radius” of the slab defined as half of its width. The distances between the inhomogeneities for these two phases are $`r_09fm`$ and $`r_213.62fm`$, respectively. Half of the oscillation period is equal to $`\frac{1}{2}\mathrm{\Delta }ϵ_F^n0.801MeV`$. On the other hand, according to the HF calculations the neutron Fermi energies are equal to:

$`ϵ_F^n(Q_p=0,\rho =0.079)`$ $`=`$ $`12.645MeV,`$ (18)
$`ϵ_F^n(Q_p=0,\rho =0.0749)`$ $`=`$ $`11.779MeV,`$ (19)
$`ϵ_F^n(Q_p=2600,\rho =0.079)`$ $`=`$ $`12.759MeV,`$ (20)
$`ϵ_F^n(Q_p=2600,\rho =0.0749)`$ $`=`$ $`11.884MeV,`$ (21)
which give for half of the periods the values:

$`{\displaystyle \frac{1}{2}}\mathrm{\Delta }ϵ_F^n(Q_p=0)`$ $`=`$ $`0.866MeV,`$ (22)
$`{\displaystyle \frac{1}{2}}\mathrm{\Delta }ϵ_F^n(Q_p=2400)`$ $`=`$ $`0.875MeV.`$ (23)
The agreement with the semiclassical estimation is very good. This demonstrates that the oscillatory effect of the energy density difference between the phases is due to the shell energy of unbound neutrons. From eq. (11), it is clear that another oscillatory effect may appear with a much longer period $`\stackrel{~}{\mathrm{\Delta }}k_F^n={\displaystyle \frac{2\pi }{\mathrm{\Delta }r_{ij}}}`$. Its mechanism is somewhat similar to the well known supershell phenomenon in atomic clusters . This effect however, requires the stability of one phase over a wide density range and therefore may be not well pronounced.

Up to now, we have neglected the pairing interaction. The pairing energy is a smooth function of the density and is not expected to play an important role for the liquid drop part of the total energy expansion. On the contrary, pairing correlations may considerably influence the shell energy .

The pairing energy can be divided with a good approximation into two parts. The first one comes from the pairing interaction acting among bound nucleons, and the second one is associated with pairing of unbound neutrons. These two contributions can be, in the first approximation, treated as independent quantities, because of the short range of the pairing interaction. For example, in the usual zero-range approximation of the pairing interaction (delta force) the amplitude for pair scattering is proportional to the integral: $`{\displaystyle \varphi _k^{}(𝐫)\varphi _l^{}(𝐫)\varphi _m(𝐫)\varphi _n(𝐫)d^3r}`$, where $`\varphi _k(𝐫)`$ are the single-particle wave functions. On the other hand, bound single-particle states, as well as resonant states inside nuclei have a small overlap with states forming continuum. Hence the integral describing the scattering of a Cooper pair from nuclei to continuum or vice versa will be relatively small, and can be, in a first approximation, neglected.

The pairing interaction acting among bound nucleons will decrease the amplitude of shell effects inside the nuclei. Note however that these effects can at most lead to a slight deformation change of nuclei and do not influence the phase transition pattern. Hence in the following we focus only on the pairing energy associated with the neutron gas. Since the pairing gap $`\mathrm{\Delta }`$ for neutron matter at subnuclear densities reaches $`34`$ MeV at the Fermi level , the occupation of the single-particle states close to the Fermi level are significantly modified by pairing. In order to analyze the influence of pairing on the shell structure, let us consider a Fermi gas of neutrons in a volume $`V`$ at a density $`\rho _n^{out}`$, interacting via a pairing force. The total energy of the system can be expressed in the BCS approximation as:

$$E=_{\mathrm{}}^{\mathrm{}}ϵv_\mu ^2(ϵ)g(ϵ)𝑑ϵ+E_{pair},$$
(25)
where $`g(ϵ)`$ is the density of single-particle states at the energy $`ϵ`$ and $`\mu `$ is the chemical potential defined by the condition:

$$\rho _n^{out}=\frac{1}{V}_{\mathrm{}}^{\mathrm{}}v_\mu ^2(ϵ)g(ϵ)𝑑ϵ.$$
(26)
The coefficient $`v_\mu ^2(ϵ)`$ represents the BCS occupation factor:

$$v_\mu ^2(ϵ)=\frac{1}{2}\left(1\frac{ϵ\mu }{\sqrt{(ϵ\mu )^2+\mathrm{\Delta }(ϵ)^2}}\right),$$
(27)
where $`\mathrm{\Delta }(ϵ)`$ is the (energy dependent) pairing gap. In order to avoid divergences, we assume that the function $`\mathrm{\Delta }`$ approaches zero sufficiently rapidly to make the integral in equation (25) convergent. The term $`E_{pair}`$ denotes the pairing energy. If the Fermi gas is subject to boundary conditions, e.g. due to presence of impurities in the system, the level density contains a fluctuating term:

$$g(ϵ,𝐝)=\stackrel{~}{g}(ϵ)+g_C(ϵ,𝐝),$$
(28)
where $`\stackrel{~}{g}(ϵ)`$ is the smooth part of the level density. The fluctuating part $`g_C(ϵ,𝐝)`$ depends on the shape and the mutual geometrical arrangement of the impurities parametrized by $`𝐝`$ and is responsible for the shell correction energy. For two spherical, cylindrical or planar obstacles located at a distance $`d`$, the semiclassical approximation gives the leading contribution to $`g_C`$ in the form (see appendix B):

$$g_C(ϵ,d)\frac{1}{4}\frac{L^iR^{2i}}{d^{1i/2}}\left(\frac{2m_n}{\pi ^2\mathrm{}^2}\right)^{\frac{i+2}{4}}ϵ^{\frac{i2}{4}}\mathrm{cos}(2kdi\frac{\pi }{4}),$$
(29)
where $`ϵ={\displaystyle \frac{\mathrm{}^2k^2}{2m_n}}`$ and $`i`$ takes the values $`0,1,2`$ for the various shapes (compare to eq. (8)).

The total energy of the system can now be divided into smooth and oscillating parts (see appendix A):

$`E`$ $`=`$ $`\stackrel{~}{E}+E_C+\mathrm{\Delta }E+E_{pair}+O(1/V)`$ (30)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}ϵv_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ)𝑑ϵ+{\displaystyle _{\mathrm{}}^{\mathrm{}}}ϵv_{\stackrel{~}{\mu }}^2(ϵ)g_C(ϵ,d)𝑑ϵ`$ (31)
$``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}v_{\stackrel{~}{\mu }}^2(ϵ)g_C(ϵ,d)𝑑ϵ{\displaystyle \frac{_{\mathrm{}}^{\mathrm{}}ϵ\frac{}{\stackrel{~}{\mu }}v_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ)𝑑ϵ}{_{\mathrm{}}^{\mathrm{}}\frac{}{\stackrel{~}{\mu }}v_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ)𝑑ϵ}}+E_{pair}+O(1/V),`$ (32)
where the first term represents the smooth part of the total energy. The term $`E_C`$ is the Casimir energy, and $`\mathrm{\Delta }E`$ provides a correction to the Casimir energy due to the Fermi energy renormalization in the presence of obstacles . The smooth part of the chemical potential $`\stackrel{~}{\mu }`$ is calculated from the condition:

$$N=_{\mathrm{}}^{\mathrm{}}v_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ)𝑑ϵ.$$
(33)
It implies that the neutron shell energy in the presence of pairing can be expressed in the form:

$$E_{shell}^{out}=E_C+\mathrm{\Delta }E+O(1/V),$$
(34)
where we have neglected corrections due to the change of pairing energy in the presence of obstacles (see appendix A).

In Fig. 5 we show the behavior of the shell energy as a function of the distance between obstacles without pairing and with a pairing interaction characterized by a pairing gap equal to $`\mathrm{\Delta }=3.5`$ MeV for $`k<20fm^1`$ (a cutoff has been introduced to avoid divergences). This pairing strength corresponds to the pairing gap at the Fermi level for a neutron Fermi momentum $`k_F^n1`$ fm-1 obtained with the Argonne $`v_{14}`$ and $`v_{18}`$ interactions . One can see that the amplitude of the shell energies are smaller but the decrease is lower than 10%, as compared to the no pairing regime. We conclude from this comparison that the pairing interaction will not change qualitatively the HF results.

## IV CONCLUSIONS

In this paper we have analyzed the structure of the neutron star crust at zero temperature using the Hartree-Fock approach with an effective Skyrme interaction. In particular, we focused on the problem of shell effects associated with unbound neutrons scattered by nuclear inhomogeneities. Our results show that:

* the relative energies of different phases fluctuate rapidly as a function of the total density,
* the fluctuations can be attributed to the shell effects associated with unbound neutrons,
* the neutron shell energy density is of the same order as the liquid drop energy density differences between the various phases present in the crust, and its behavior can be easily understood in terms of the periodic orbits in the system,
* the pairing correlations lead only to a slight decrease of the shell energy and are not expected to influence significantly the structure of the crust.

Consequently, we may conclude that the structure of the crust may be quite complicated since the shell effects associated with unbound neutrons may easily reverse the phase transition order predicted by liquid drop based approaches. Moreover the number of phase transitions may increase since the same phase may appear for various density ranges (see Fig. 4).

Our results suggest also that several phases, differing by the lattice type, could coexist in the crust. This possibility was not taken into account in previous investigations since in a liquid drop based approach, the system favors only one lattice type for a given nuclear shape. However, since the neutron shell energy is very sensitive to the spatial order in the system, this simple picture will change, and one may expect the coexistence of various lattice geometries with different lattice constants. Note also that once a phase is formed, its geometric order will be stabilized by both the Coulomb repulsion and the shell energy. Indeed, although the Coulomb energy is a smooth function of the nuclear displacement, the shell energy is not. Since several different orbits contribute to the shell effects (except for the slab–like phase) the displacement of a single nucleus from its equilibrium position in the lattice will give rise to interference effects depending on the lattice type. Hence it is likely that the system will favor distorted lattices, or lattices with defects which decrease the shell energy. In that respect the present results support the conclusions drawn on the basis of a simple model .

All these arguments suggest that the inner crust is an extremely rich and complicated system which may turn out to be on the verge of a disordered phase.

###### Acknowledgements.

This research was supported in part by the Polish Committee for Scientific Research (KBN) under Contract No. 5 P03B 014 21 and the Wallonie/Brussels-Poland integrated action program. Numerical calculations were performed at the Interdisciplinary Centre for Mathematical and Computational Modelling (ICM) at Warsaw University. Authors would like to thank Hubert Flocard for providing us with the numerical code solving the Poisson equation and Aurel Bulgac for discussions and valuable comments.

## A Shell energy of the Fermi gas with impurities

We derive here expressions for the shell energy associated with impurities immersed in a fermionic environment.

Let us consider a Fermi gas confined in a volume $`V`$, large compared to both the Fermi wavelength ($`k_F^3V>>1`$) and the volume of the obstacles. At the end of the calculations the limit $`V\mathrm{}`$ will be taken. The shell energy is an interaction energy between “foreign” objects immersed in the gas. Its origin comes from the modification of the level density by the presence of obstacles. The total level density can be decomposed in the following way :

$$g(ϵ,𝐝)=\stackrel{~}{g}(ϵ)+g_C(ϵ,𝐝),$$
(A1)
where $`\stackrel{~}{g}`$ is the smooth part of $`g`$ and can be obtained either from the Weyl expansion  (in the case of impenetrable objects) or within the Thomas-Fermi approach . The set of parameters $`𝐝=(d_1,d_2,\mathrm{},d_n)`$ describes the geometrical arrangements of the obstacles. The smooth part can also be treated as the total density of states for a system with obstacles located at infinity. Namely:

$$\underset{d_i\mathrm{}}{lim}g(ϵ,𝐝)=\stackrel{~}{g}(ϵ),\text{ for i=1,2,…n}.$$
(A2)
The fluctuating part of the level density $`g_C`$ is associated with the formation of resonances. Its evaluation in a semiclassical approximation is presented in the next section.

The shell energy can be defined in a form similar to the Casimir energy :

$$E_{shell}(𝐝)=_{\mathrm{}}^\mu 𝑑ϵϵg(ϵ,𝐝)_{\mathrm{}}^{\stackrel{~}{\mu }}𝑑ϵϵ\stackrel{~}{g}(ϵ),$$
(A3)
where $`\mu `$ and $`\stackrel{~}{\mu }`$ denote the chemical potentials associated with the total and smooth level densities, respectively. They are related to each other through the particle number conservation condition:

$$N=_{\mathrm{}}^\mu 𝑑ϵg(ϵ,𝐝)=_{\mathrm{}}^{\stackrel{~}{\mu }}𝑑ϵ\stackrel{~}{g}(ϵ).$$
(A4)
The decomposition of the integrals in equation (A3) results in the following expression:

$$E_{shell}(𝐝)=\mathrm{\Delta }E+E_C=_{\stackrel{~}{\mu }}^\mu ϵ\stackrel{~}{g}(ϵ)𝑑ϵ+_{\mathrm{}}^\mu ϵg_C(ϵ,𝐝)𝑑ϵ,$$
(A5)
where $`\mathrm{\Delta }E`$ represents the correction to the shell energy due to the variation of the chemical potential induced by the presence of obstacles, while $`E_C`$ corresponds to the ”pure Casimir–like” force. Since $`\mathrm{\Delta }\mu =\mu \stackrel{~}{\mu }O(1/V)`$ and $`\stackrel{~}{g}V`$ then

$`{\displaystyle _{\stackrel{~}{\mu }}^\mu }ϵ\stackrel{~}{g}(ϵ)𝑑ϵ`$ $`=`$ $`\mu \mathrm{\Delta }\mu \stackrel{~}{g}(\mu )+O(1/V),`$ (A6)
$`{\displaystyle _{\stackrel{~}{\mu }}^\mu }\stackrel{~}{g}(ϵ)𝑑ϵ`$ $`=`$ $`\mathrm{\Delta }\mu \stackrel{~}{g}(\mu )+O(1/V),`$ (A7)
and thus it is clear that both $`\mathrm{\Delta }E`$ and $`E_C`$ do not vanish in the limit of infinite volume. The condition (A4) can be used to calculate $`\mathrm{\Delta }\mu `$:

$$\mathrm{\Delta }\mu =\frac{1}{\stackrel{~}{g}(\mu )}_0^\mu g_C(ϵ,𝐝)𝑑ϵ+O(1/V^2).$$
(A8)
The shell energy can be finally written in the form:

$$E_{shell}(𝐝)=_{\mathrm{}}^\mu (ϵ\mu )g_C(ϵ,𝐝)𝑑ϵ+O(1/V).$$
(A9)
Taking the limit $`V\mathrm{}`$, the higher order terms vanish and the first term in the above formula yields an exact expression for the shell energy.

In the case of a superfluid Fermi gas, the energy gap depends on the momentum $`k`$. To avoid divergences in the expression for the total energy, one has to assume that $`\mathrm{\Delta }(k)`$ is tending to zero at least like $`1/k`$ for $`k\mathrm{}`$. The pairing properties of the system are affected by the presence of obstacles and thus the solution of the gap equation with impurities $`\mathrm{\Delta }(k,𝐝)`$ differs slightly from the smooth pairing gap $`\stackrel{~}{\mathrm{\Delta }}(k)`$ obtained using the smooth level density: $`\mathrm{\Delta }(k,𝐝)=\stackrel{~}{\mathrm{\Delta }}(k)+\mathrm{\Delta }_C(k,𝐝)`$. $`\mathrm{\Delta }_C`$ denotes the quantum correction to the smooth pairing gap and is due to the fact that the level density is modified. Hence, eq. (A3) is generalized for superfluid systems in the form:

$`E_{shell}(𝐝)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ϵϵv_\mu ^2(ϵ,𝐝)g(ϵ,𝐝){\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ϵϵ\stackrel{~}{v}_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ)`$ (A10)
$`+`$ $`E_{pair}(𝐝)\stackrel{~}{E}_{pair},`$ (A11)
where $`v_\mu ^2`$ is the BCS occupation factor in the presence of obstacles (27) and the term $`\stackrel{~}{v}_{\stackrel{~}{\mu }}`$ represents the smooth occupation factor obtained for impurities far apart from each other. The last two terms are due to the correction of the pairing energy related to the appearance of resonant states between obstacles.

For simplicity, we assume that $`\mathrm{\Delta }_C=0`$ (and thus $`\stackrel{~}{v}_\mu ^2=v_\mu ^2`$) and neglect the effect of the geometrical arrangements of the obstacles on the pairing properties of the system. The equation simplifies then to the form:

$$E_{shell}(𝐝)=_{\mathrm{}}^{\mathrm{}}𝑑ϵϵv_\mu ^2(ϵ)g(ϵ,𝐝)_{\mathrm{}}^{\mathrm{}}𝑑ϵϵv_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ),$$
(A12)
where the chemical potentials $`\mu `$ and $`\stackrel{~}{\mu }`$ are determined from the conditions:

$$N=_{\mathrm{}}^{\mathrm{}}𝑑ϵv_\mu ^2(ϵ)g(ϵ,𝐝)=_{\mathrm{}}^{\mathrm{}}𝑑ϵv_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ).$$
(A13)
Expanding $`v_\mu ^2`$ around $`\stackrel{~}{\mu }`$ one gets:

$$E_{shell}(𝐝)=\mathrm{\Delta }\mu _{\mathrm{}}^{\mathrm{}}𝑑ϵϵ\frac{v_{\stackrel{~}{\mu }}^2(ϵ)}{\stackrel{~}{\mu }}\stackrel{~}{g}(ϵ)+_{\mathrm{}}^{\mathrm{}}𝑑ϵϵ\left(v_{\stackrel{~}{\mu }}^2(ϵ)+\mathrm{\Delta }\mu \frac{v_{\stackrel{~}{\mu }}^2(ϵ)}{\stackrel{~}{\mu }}\right)g_C(ϵ,𝐝)+O(1/V).$$
(A14)
Note however that the term: $`\mathrm{\Delta }\mu {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ϵϵ{\displaystyle \frac{v_{\stackrel{~}{\mu }}^2(ϵ)}{\stackrel{~}{\mu }}}g_C(ϵ,𝐝)`$ is also proportional to $`1/V`$ and thus vanishes when $`V\mathrm{}`$. Determining $`\mathrm{\Delta }\mu `$ from condition (A13) one gets the expression for the shell energy:

$$E_{shell}(𝐝)=_{\mathrm{}}^{\mathrm{}}ϵv_{\stackrel{~}{\mu }}^2(ϵ)g_C(ϵ,𝐝)𝑑ϵ_{\mathrm{}}^{\mathrm{}}v_{\stackrel{~}{\mu }}^2(ϵ)g_C(ϵ,𝐝)𝑑ϵ\frac{_{\mathrm{}}^{\mathrm{}}ϵ\frac{}{\stackrel{~}{\mu }}v_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ)𝑑ϵ}{_{\mathrm{}}^{\mathrm{}}\frac{}{\stackrel{~}{\mu }}v_{\stackrel{~}{\mu }}^2(ϵ)\stackrel{~}{g}(ϵ)𝑑ϵ}+O(1/V).$$
(A15)
Since $`v_{\stackrel{~}{\mu }}^2(ϵ)\theta (\stackrel{~}{\mu }ϵ)`$ and $`{\displaystyle \frac{v_{\stackrel{~}{\mu }}^2(ϵ)}{\stackrel{~}{\mu }}}\delta (ϵ\stackrel{~}{\mu })`$ when $`\mathrm{\Delta }(k)`$ tends to zero, the equation (A15) becomes identical with (A3) in the no pairing limit.

## B Evaluation of the fluctuating part of level density and the shell energy in the semiclassical approximation

We derive here expressions for the shell (interaction) energy for two identical, impenetrable obstacles at large separations (8). Three geometries of the obstacles are considered: spherical, rod–like and planar geometries. The system is characterized by a set of parameters $`R,L,d`$ which denote the radius, length (in the case of rods and slabs), and the distance between obstacles, respectively. The distance is measured as half of the length of the shortest periodic orbit. It means that e.g. in the case of spheres the distance between their centers is equal to $`d+2R`$. In the case of a slab, $`R`$ denotes half of its width.

In the semiclassical approximation, the fluctuating part of the level density $`g_C`$ can be obtained through the summation of the contributions associated with the classically allowed periodic orbits in the system. Recent exact calculations based on the S-matrix approach indicated that this approximation work surprisingly well even if the obstacles are close to each other, i.e. when $`d<<R`$. In general, the summation is a tedious task if many orbits contribute to the level density. Especially when elliptic or parabolic orbits are considered, the problem requires a more careful treatment in order to avoid divergences . The semiclassical approximation to the interaction energy between two obstacles is particularly useful if the periodic orbits of the system are unstable (hyperbolic) and isolated . It holds true for spherical obstacles. For two spheres of radius $`R`$, located at a distance $`d`$, the contribution to the level density arising from the single periodic orbit is:

$$g_C^{sphere}(k,d)=\frac{md}{\pi \mathrm{}^2k}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{2\mathrm{sinh}^2(\frac{n\lambda }{2})}\mathrm{cos}(2nkd\sigma \frac{\pi }{2}),$$
(B1)
where $`\lambda =2\mathrm{log}\left(1+{\displaystyle \frac{d}{R}}+\sqrt{{\displaystyle \frac{d}{R}}\left({\displaystyle \frac{d}{R}}+2\right)}\right)`$ is the Lyapunov exponent determining the stability of the orbit, $`\sigma `$ is the Maslov index depending on the type of boundary conditions, and $`ϵ={\displaystyle \frac{\mathrm{}^2k^2}{2m}}`$. For Dirichlet boundary conditions, which we assume here, each reflection changes the phase by $`\pi `$, thus in our case $`\sigma =4n`$ where $`n`$ is the number of repetitions of periodic orbit. The shell energy calculated from the expression (A3) gives :

$$E_{shell}^{sphere}(d)=\frac{\mathrm{}^2}{16\pi md^2}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^3\mathrm{sinh}^2(\frac{n\lambda }{2})}\left(2nk_Fd\mathrm{cos}(2nk_Fd)\mathrm{sin}(2nk_Fd)\right).$$
(B2)
At large distances, i.e. for $`{\displaystyle \frac{d}{R}}>>1`$, the sum can be truncated at $`n=1`$ and the expression simplifies to:

$$E_{shell}^{sphere}(d)\frac{\mathrm{}^2R^2}{8m}\left(\frac{3}{\pi }\right)^{1/3}\frac{\rho ^{1/3}}{d^3}\mathrm{cos}(2k_Fd)+O(1/d^5),$$
(B3)
where we have expressed the amplitude of the shell energy in terms of the density $`\rho ={\displaystyle \frac{k_F^3}{3\pi ^2}}`$. This expression is associated with the fluctuating level density:

$$g_C^{sphere}(k,d)\frac{1}{2}\frac{mR^2}{\pi \mathrm{}^2kd}\mathrm{cos}(2kd)+O(1/d^3).$$
(B4)
In the case of two parallel rods of length $`L`$, the orbit is no longer isolated. Hence it is convenient to consider first the two-dimensional case of two disks. Then the single periodic orbit gives rise to the fluctuating part of the level density:

$$g_C^{disk}(k,d)=\frac{md}{\pi \mathrm{}^2k}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{\mathrm{sinh}(\frac{n\lambda }{2})}\mathrm{cos}(2nkd).$$
(B5)
In order to obtain the $`g_C^{rod}`$ one has to calculate first the fluctuating part of the particle number $`N_C^{rod}(\mu ,d)={\displaystyle _{\mathrm{}}^\mu }g_C^{rod}(ϵ,d)𝑑ϵ`$. But $`N_C^{rod}`$ has to be an integral of the product of $`g_C^{disk}`$ and the level density of one-dimensional Fermi gas. Namely, the quantity $`N_C^{rod}`$ can be expressed as:

$`N_C^{rod}(k_F,d)`$ $`=`$ $`{\displaystyle \frac{L}{2\pi }}{\displaystyle _0^{k_F}}g_C^{disk}(k)d\left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}\right){\displaystyle _{\sqrt{k_F^2k^2}}^{\sqrt{k_F^2k^2}}}𝑑k_p`$ (B6)
$`=`$ $`{\displaystyle \frac{Lk_F}{4\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n\mathrm{sinh}(\frac{n\lambda }{2})}}J_1(2nk_Fd),`$ (B7)
where $`k_p`$ denotes the momenta parallel to rods and $`J_1`$ is the Bessel function of the first kind. Since $`g_C^{rod}(k,d)={\displaystyle \frac{dN_C^{rod}}{dk}}{\displaystyle \frac{m}{\mathrm{}^2k}}`$ we get***In the Refs.  the factor “$`1/\pi `$” has been missed in the expression for the shell energy and level density of the rod–like phase.:

$$g_C^{rod}(k,d)=\frac{mL}{4\pi \mathrm{}^2k}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n\mathrm{sinh}(\frac{n\lambda }{2})}\left(kdJ_0(2nkd)+\frac{1}{n}J_1(2nkd)kdJ_2(2nkd)\right),$$
(B8)
and the shell energy:

$$E_{shell}^{rod}(d)=\frac{\mathrm{}^2k_F^2}{8m\pi }\frac{L}{d}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2\mathrm{sinh}(\frac{n\lambda }{2})}J_2(2nk_Fd).$$
(B9)
In the limit $`{\displaystyle \frac{d}{R}}>>1`$, using the asymptotic formula for a Bessel function and expressing the Fermi momentum as a function of the density, one gets:

$$g_C^{rod}(k,d)\frac{1}{2}\frac{mLR}{\sqrt{kd\pi ^3}\mathrm{}^2}\mathrm{cos}(2kd\frac{\pi }{4})+O(1/d^{3/2}),$$
(B10)
and for the shell energy:

$$E_{shell}^{rod}(d)\frac{\mathrm{}^2LR}{8m}\left(\frac{3}{\pi }\right)^{1/2}\frac{\rho ^{1/2}}{d^{5/2}}\mathrm{cos}(2nk_Fd\frac{\pi }{4})+O(1/d^{7/2}).$$
(B11)
In the case of two slabs the semiclassical approximation is no longer applicable since the classical orbits in this case are parabolic which implies that $`\lambda =0`$ and the expression for $`g_C`$ diverges. On the other hand, the problem can be solved analytically yielding an exact solution for the shell energy . In the following, however we derive an expression in a different way, which has the advantage that the final formula has a form similar to the semiclassical expressions for spheres and rods.

We consider first the fluctuating level density of a one-dimensional Fermi gas confined in a box of length $`d`$ with Dirichlet boundary conditions. The total level density in this case can be calculated exactly and has the form :

$$g(k,d)=\frac{md}{\mathrm{}^2k\pi }(1+2\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{cos}(2nkd)).$$
(B12)
The second term corresponds to the fluctuating part of the level density:

$$g_C^{1dim}(k,d)=\frac{2md}{\mathrm{}^2k\pi }\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{cos}(2nkd).$$
(B13)
Using an analogous approach as in the case of rods, one can calculate the fluctuating part of the particle number as a function of the momentum:

$`N_C^{slab}(k_F)`$ $`=`$ $`\left({\displaystyle \frac{L}{2\pi }}\right)^2{\displaystyle _0^{k_F}}g_C^{1dim}(k_z,d)d\left({\displaystyle \frac{\mathrm{}^2k_z^2}{2m}}\right){\displaystyle _{k_x^2+k_y^2<k_F^2k_z^2}}𝑑k_x𝑑k_y`$ (B14)
$`=`$ $`\left({\displaystyle \frac{L}{2\pi }}\right)^2{\displaystyle \frac{1}{2d^2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2nk_Fd\mathrm{cos}(2nk_Fd)+\mathrm{sin}(2nk_Fd)}{n^3}},`$ (B15)
where $`k_z`$ is the momentum component perpendicular to slabs. Hence the fluctuating level density for slabs is given by:

$$g_C^{slab}(k,d)=\frac{dN_C^{slab}}{dk}\frac{m}{\mathrm{}^2k}=2m\left(\frac{L}{2\pi \mathrm{}}\right)^2\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}(2nkd)}{n}.$$
(B16)
The shell energy calculated using expression (A3) gives:

$$E_{shell}^{slab}=\frac{\mathrm{}^2L^2d}{32m\pi ^2}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{(nd)^5}(6nk_Fd\mathrm{cos}(2nk_Fd)+4(nk_Fd)^2\mathrm{sin}(2nk_Fd)3\mathrm{sin}(2nk_Fd)).$$
(B17)
Note that unlike in the previous cases the truncation of the summation for large separation is not justified since all terms have the same dependence on the distance. Nevertheless, the gross shell structure is still represented by the $`n=1`$-term and thus for comparison with previous expressions we consider the leading term $`n=1`$:

$$E_{shell}^{slab}(d)\frac{\mathrm{}^2L^2}{8m}\left(\frac{9}{\pi ^2}\right)^{1/3}\frac{\rho ^{2/3}}{d^2}\mathrm{sin}(2k_Fd)+O(1/d^2),$$
(B18)
and for the fluctuating level density:

$$g_C^{slab}(k,d)\frac{1}{2}\frac{mL^2}{\pi ^2\mathrm{}^2}\mathrm{cos}(2kd\frac{\pi }{2})+O(1/d^0).$$
(B19)
The expressions (B3,B11, B18) and (B4,B10,B19) can be gathered into a single formula for the shell energy and the fluctuating part of the level density, respectively:

$$E_{shell}^i(d)\frac{\mathrm{}^2L^iR^{2i}}{8m}\left(\frac{3}{\pi }\right)^{\frac{2+i}{6}}\frac{(\rho )^{\frac{2+i}{6}}}{d^{\frac{6i}{2}}}\mathrm{cos}\left(2k_Fdi\frac{\pi }{4}\right)+O(1/d^{5\frac{3i}{2}}),$$
(B20)
$$g_C^i(ϵ,d)\frac{1}{4}\frac{L^iR^{2i}}{d^{1i/2}}\left(\frac{2m}{\pi ^2\mathrm{}^2}\right)^{\frac{i+2}{4}}ϵ^{\frac{i2}{4}}\mathrm{cos}\left(2kdi\frac{\pi }{4}\right)+O(1/d^{3\frac{3i}{2}}),$$
(B21)
where $`i=0,1,2`$ correspond to the spherical, rod and slab case, respectively and $`ϵ={\displaystyle \frac{\mathrm{}^2k^2}{2m}}`$.


format: text

text:
# TUM-HEP-449/01 Effects of random matter density fluctuations on the neutrino oscillation transition probabilities in the Earth

## 1 Introduction

The effects of matter on neutrino oscillations  in the Earth have been investigated in various contexts and with several models . It is now well-known that the matter density can significantly change the reconstructed neutrino energy spectrum produced by long baseline neutrino experiments, such as by neutrino factories . For most calculations the Preliminary Reference Earth Model (PREM) density profile  has been used, which is obtained from geophysical seismic wave measurements (see, e.g., Refs.  for information about the structure of the Earth’s interior). Furthermore, small errors in the PREM matter density with up to 5% amplitude have been found and documented by many geophysics groups (for a summary, see, e.g., Ref. ). Note, however, that the Earth’s matter density distribution is not directly observable from seismological data . In this paper, we will discuss how these fluctuations in the Earth matter density affect the neutrino oscillation transition probabilities.

It has been found that for short baselines matter effects are small . Thus, any fluctuations could be treated as second order effects and can therefore be neglected. For long baselines, however, the fluctuations may be significant. It was noted in Ref.  that this effect can be important as an additional uncertainty in the determination of the CP phase $`\delta _{CP}`$, especially for certain values of $`\delta _{CP}`$. To estimate this effect, the authors used a logarithmic distribution with a certain length scale and amplitude with a path integral method for the numerical evaluation. On the other hand, it was shown in Ref. , using a perturbation theoretical approach, that fluctuations with small amplitudes on length scales much shorter than the oscillation length in matter average out and give no net effect at all. In this paper, we are interested in the errors on the transition probabilities as functions of the length scale and amplitude of the matter density fluctuations. We will, in particular, focus on the errors on the appearance probability of electron neutrinos and electron antineutrinos at typical neutrino factory energies, since matter density noise effects could be rather substantial in the determination of the CP phase and matter effects are largest in the appearance channel.

## 2 The effect of a matter density perturbation

Before we come to modeling the fluctuations in the Earth matter density profile, let us study the effect of a single perturbation in the matter density. In any quantum mechanical system, described by a Schrödinger equation, the impact on a free particle’s motion of a potential depends on its length scale as well as its amplitude. When the length scale of the potential is much shorter than the characteristic wave length of the incident particle, the particle’s wave function will be unable to resolve the exact spatial structure of this potential. It is then possible to replace the potential by a $`\delta `$-distribution with an amplitude equal to the integral of the original potential. Similarly, in neutrino oscillations, described by a Schrödinger equation, a perturbation on a length scale much shorter than the oscillation length in matter could be replaced by a $`\delta `$-distribution.

This can be seen in Fig. 1 (left plot), where the relative error $`|\mathrm{\Delta }P_{\mu e}|/P_{\mu e}^{\mathrm{ref}}|P_{\mu e}P_{\mu e}^{\mathrm{ref}}|/P_{\mu e}^{\mathrm{ref}}`$, coming from a rectangular matter density perturbation, is plotted as a function of the product of the length scale $`\lambda /L_{\mathrm{matter}}^{\mathrm{osc}}`$ and the density contrast $`\mathrm{\Delta }\rho /\overline{\rho }`$. For the oscillation parameters we choose $`\theta _{12}=45^{}`$, $`\theta _{23}=45^{}`$, $`\theta _{13}=5^{}`$, $`\mathrm{\Delta }m_{32}^2=2.510^3\mathrm{eV}^2`$, $`\mathrm{\Delta }m_{21}^2=3.6510^5\mathrm{eV}^2`$, and $`\delta _{CP}=0`$, corresponding to the LMA solution with a value for $`\theta _{13}`$ somewhat below the CHOOZ bound neglecting CP violating effects. The other parameters (and terms) are given and described in the figure caption. The relative error is shown for several fixed values of $`\lambda L_{\mathrm{matter}}^{\mathrm{osc}}`$, i.e., only the matter density is varied in such a way that the area of the perturbation is unaffected for a fixed value on the horizontal axis. For the cut at the vertical line this constant area of the matter density perturbation is shown in the right plot, indicating that the corresponding parameter value is already a quite pessimistic choice. Apparently, the curves in the left plot are approximately equal to each other at least below the parameter value indicated by the vertical line and especially for very small fixed values of $`\lambda `$. This means that the relevant parameter for $`\lambda L_{\mathrm{matter}}^{\mathrm{osc}}`$ is the integral of the matter density perturbation, i.e., the product of the length scale and the amplitude. Thus, we may expect that interference effects, arising from non-commuting operators in the Hamiltonian corresponding to different matter density layers, become irrelevant for very short length scales compared with the oscillation length in matter. This can also be seen in the analytical perturbation theoretical approach in App. A, where the case of baselines much shorter than the oscillation length leads to first order corrections $`\lambda \mathrm{\Delta }\rho `$ to the transition probabilities, i.e., the product of the length scale and the amplitude is the relevant parameter. Another result of the numerical analysis is that we do not have to take into account isolated short scale perturbations, such as the ones coming from the matter density contrast in, for example, a mine. Estimating the length scale of such a perturbation to be shorter than $`10\mathrm{km}`$ and the relative density contrast to be of the order of $`10\%`$, we can read off a relative error much smaller than $`1\%`$ from Fig. 1. Since the amplitude of the appearance probability $`P_{\mu e}`$ is basically proportional to $`\mathrm{sin}^22\theta _{13}`$ and the error in the determination of $`\theta _{13}`$ is rather substantial , we will henceforth neglect very short isolated perturbations.

## 3 A model for the matter density fluctuations

In this section, we will construct a model for the fluctuations in the Earth matter density. Figure 2 shows the percentage fluctuations in the Earth matter density at a depth of $`20\mathrm{km}`$ below the Earth’s surface obtained from seismic wave measurements.

Though a neutrino beam traverses a large range of different depths, one may take this figure as an estimate of the characteristic length scales and amplitudes involved in the problem. It suggests length scales of the order of some thousands of kilometers and amplitudes of the order of $`\pm 4\%`$, whereas at greater depths one can indeed have somewhat larger amplitudes. In addition, the length scales and amplitudes do not seem to vary too much around their average absolute values. This also implies that the transition regions between negative and positive amplitudes are quite short compared to the overall structure. So why not simply use these measurements in neutrino physics, instead of discussing uncertainties in the Earth matter density? First, these measurements contain some averaging as well as uncertainties in the equation of state of the Earth matter density profile from the seismic wave velocity profile (see, e.g., Refs. ). Second, different groups obtain different results , which are, however, not qualitatively so much different with respect to the parameters we will identify below.

In order to model these fluctuations realistically and investigate the dependence of the relevant parameters on the neutrino oscillation transition probabilities, we use a step-function approach, varying the absolute value of the amplitude $`\mathrm{\Delta }\rho >0`$ and the length scale $`\lambda >0`$ around some average values $`\mathrm{\Delta }\rho _0>0`$ and $`\lambda _0>0`$ at random. For the random variation we choose a Gaussian distribution with standard deviations $`\sigma _\lambda `$ and $`\sigma _{\mathrm{\Delta }\rho }`$, respectively, which is truncated at zero. Figure 3 shows some sample profiles for different values of these standard deviations and for fixed average values.

Comparing this model with Fig. 2, we can then estimate a set of realistic parameters as

$`\lambda _0`$ $``$ $`2000\mathrm{km},`$ (1a)
$`\sigma _\lambda `$ $``$ $`1500\mathrm{km}=0.75\lambda _0,`$ (1b)
$`\mathrm{\Delta }\rho _0`$ $``$ $`3\%\overline{\rho },`$ (1c)
$`\sigma _{\mathrm{\Delta }\rho }`$ $``$ $`1\%\overline{\rho }=1/3\mathrm{\Delta }\rho _0,`$ (1d)
where $`\overline{\rho }=3\mathrm{g}/\mathrm{cm}^3`$ is the average matter density.

## 4 Numerical analysis

For physical reasons we do not expect surprises in certain regions of the parameter space. Since we cannot show the results for the whole parameter space in $`\lambda _0`$, $`\mathrm{\Delta }\rho _0`$, $`\sigma _\lambda `$, and $`\sigma _{\mathrm{\Delta }\rho }`$ simultaneously, we will now systematically investigate the dependence on some parameters by keeping the other ones fixed. One could do this by showing either the absolute errors in the appearance probabilities coming from matter density fluctuations $`|\mathrm{\Delta }P_{\alpha \beta }|`$, or the relative errors $`|\mathrm{\Delta }P_{\alpha \beta }|/P_{\alpha \beta }^{\mathrm{ref}}`$. For the neutrino channel $`\nu _\mu \nu _e`$, the transition probabilities are in most regions relatively large. In this case, the relative errors are quite meaningful and are usually some percent of the total probabilities. Depending on the parameters they can sometimes even be larger than $`10\%`$. However, for the antineutrino channel $`\overline{\nu }_\mu \overline{\nu }_e`$ the absolute probabilities in the denominators of the relative errors are rather small. Therefore, it turns out that the relative errors are not very sensible in this case. Comparing plots for the absolute and relative errors and mainly focusing on the qualitative behavior, we thus decided to show only the absolute error plots. Nevertheless, the fact that the absolute values of these errors are rather small does not mean that they are small compared to the transition probabilities. For the simulations, a large number of matter density profiles was created at random and the relative error was averaged over all computations with these profiles.

Figure 4 shows the absolute errors in the appearance probabilities for the neutrino channel, $`P_{\mu e}`$, and the antineutrino channel, $`P_{\overline{\mu }\overline{e}}`$, plotted as functions of $`\sigma _\lambda `$ for a fixed value of $`\sigma _{\mathrm{\Delta }\rho }`$ and $`\sigma _{\mathrm{\Delta }\rho }`$ for a fixed value of $`\sigma _\lambda `$, respectively. For the fixed parameter values we choose, if not otherwise noted, the values from Eq. (1).

In all of our plots, the vertical lines correspond to the parameter values in Eq. (1), i.e., the same point in the multi-dimensional parameter space. Since the curves in Fig. 4 are slowly varying close to the vertical lines, we will further on take the respective values from Eq. (1) for the standard deviations. In general, the absolute errors are growing with larger fluctuations in the length scale and the amplitude. This means, in our model, that for more irregular matter density fluctuations we obtain larger errors in the transition probabilities. A possible explanation could be the zero-truncations of the Gaussian distributions, which means that for large standard deviations the average values are shifted. For antineutrinos the absolute errors are in all cases much smaller than for neutrinos, because for antineutrinos matter effects are, in general, much smaller for energies larger than a few GeV (no resonance effects).

Next, let us investigate the dependence of the absolute errors on the length scale $`\lambda _0`$ and the amplitude $`\mathrm{\Delta }\rho _0`$. Here we choose the relative standard deviations to be $`\sigma _\lambda =0.75\lambda _0`$ and $`\sigma _{\mathrm{\Delta }\rho }=1/3\mathrm{\Delta }\rho _0`$, corresponding to Eq. (1) for the appropriate $`\lambda _0`$ and $`\mathrm{\Delta }\rho _0`$. The result of this analysis is shown in Fig. 5, from which it can be seen that the absolute error is essentially proportional to the length scale as well as the amplitude of the fluctuations.

Indeed, this result is, with respect to the average length scale and amplitude, very similar to what we obtained from a single perturbation in the Earth matter density. In the latter case, we observed that the product of the length scale and amplitude determines the error in the probabilities. Thus, fixing one of these two parameters gives a linear dependence on the other one. This is, of course, only true for a single perturbation without interference effects, which means that corrections to linearity have to be taken into account in our more general model. One such interesting correction is the bumps in the upper plot. Since in our model the number of steps $`n`$, in which the matter density profile has been divided into, depends on the length scale, for a small number of steps $`n`$ the (average) transition $`nn1`$ can be seen as a bump in the absolute errors. For large $`n`$ the relative contribution of this effect becomes negligible.

Finally, taking the parameter values as given in Eq. (1), the absolute errors in the appearance probabilities $`P_{\mu e}`$ and $`P_{\overline{\mu }\overline{e}}`$ are plotted as functions of energy and baseline in Fig. 6.

From the upper plot we can observe that the absolute errors are rather small for large energies ($`E10\mathrm{GeV}`$) and so are the relative errors at least for neutrinos. Since for a neutrino factory the energy spectrum for muon neutrinos is peaked at the maximum energy , the low-energy contributions contain less statistical information. This point supports our focus on high energies in most of the plots.

As far as the baseline dependence in the lower plot of Fig. 6 is concerned, the absolute errors are vanishing for small baselines, because matter effects become, in general, negligible. Although the neutrino and antineutrino curves look quite different in this plot, the basic principle of the shapes is the same. They are both periodic functions of the baseline, as can easily be seen for the antineutrinos, and the difference is only due to different period and amplitude.

## 5 Summary and conclusions

In this paper, we first investigated the effects of a single perturbation in the average matter density of the Earth (especially the Earth’s mantle). We then introduced a model for the description of random matter density fluctuations in the Earth’s mantle, which is based on the observations obtained from geophysical measurements. Finally, we used this model to analyze the dependence of the absolute errors in the appearance probabilities $`P_{\mu e}`$ and $`P_{\overline{\mu }\overline{e}}`$, especially important for long baseline neutrino factory CP measurements, on our model parameters.

We observed that, in particular for neutrinos, the absolute error in the appearance probability for random matter density fluctuations can be quite substantial, corresponding to relative errors of some few percent. The error is essentially directly proportional to the product of the amplitude and length scale of a single perturbation or the random fluctuations. Furthermore, for short baselines $`L1000\mathrm{km}`$, the errors are vanishing for both neutrinos and antineutrinos which means that matter density fluctuations can in this case be regarded as second order effects compared to the small matter effects.

Finally, we comment on the three possible scenarios suggested in the summary of Ref. :

1. The uncertainty of present density models poses no significant problems.
2. Moderate reduction of the uncertainty, through more detailed analysis of data, is required.
3. Significant reduction of this uncertainty, by conducting a large scale campaign of geophysical observations, is required.

From our analysis, scenario (2) best fits our conclusions. In our calculations with randomly generated matter density profiles, we assumed that the fluctuations are completely unknown and we obtained relative errors of the order of magnitude of some few percent. Comparing the results of the measurements of the PREM profile corrections of different groups  indicates that there is not yet sufficient agreement on the data. Moderate reduction of the uncertainty by a more detailed analysis of data should help to settle this problem for neutrino physics.

## Appendix A A perturbation theoretical approach to a single perturbation in the Earth matter density

In perturbation theory, we can show that for a single perturbation with constant amplitude $`\mathrm{\Delta }A\pm 2\sqrt{2}G_FE\mathrm{\Delta }n_e\mathrm{\Delta }\rho `$ in the average Earth matter density $`A\overline{\rho }`$, the first order perturbation term is proportional to the “area” $`\mathrm{\Delta }S\lambda \mathrm{\Delta }A\lambda \mathrm{\Delta }\rho `$ of the perturbation, i.e., the length scale $`\lambda `$ times the amplitude $`\mathrm{\Delta }A`$.

The Hamiltonian for the propagation of the neutrinos is

$$(r)=_0+A𝒦+\mathrm{\Delta }A𝒦+\mathrm{\Delta }A𝒦,$$
(2)
where $`𝒦|\nu _e\nu _e|`$ is the projector onto the flavor state $`|\nu _e`$. Let $`𝒱(r)=e^{ir}=_{a=1}^3e^{i\xi _ar}𝒫_a`$ be the unperturbed evolution operator, where $`𝒫_a|\nu _a\nu _a|`$ is the projector onto the mass eigenstate $`|\nu _a`$ with the eigenvalue $`\xi _a`$ in matter. Then the evolution operator of the full evolution equation can, to first order in perturbation theory, be written as

$$𝒰(r)𝒱(r)i𝒱(r)\mathrm{\Delta }A\underset{a=1}{\overset{3}{}}\underset{b=1}{\overset{3}{}}f_{ab}𝒫_a𝒦𝒫_b,$$
(3)
where

$$f_{ab}\lambda e^{i(\xi _a\xi _b)r_0}\frac{\mathrm{sin}\left(\xi _a\xi _b\right)\frac{\lambda }{2}}{\left(\xi _a\xi _b\right)\frac{\lambda }{2}}$$
(4)
and $`r_0`$ is the position of the center of the fluctuation. For a fluctuation length $`\lambda `$, which is much shorter than the oscillation length, i.e., $`(\xi _a\xi _b)(\lambda /2)1`$, the last factor in Eq. (4) is approximately equal to unity. We then obtain

$$f_{ab}\lambda e^{i(\xi _a\xi _b)r_0}.$$
(5)
With this result it is easy to see that the transition probability $`P_{\alpha \beta }(r)`$ at the position $`r`$ can be written as

$$P_{\alpha \beta }(r)P_{\alpha \beta }^0(r)+2\mathrm{\Delta }S\mathrm{}[X^{}(r)Y(r)],$$
(6)
where $`P_{\alpha \beta }^0(r)`$ is the unperturbed transition probability in constant matter density, $`X(r)\nu _\beta |𝒱(r)|\nu _\alpha `$, and $`Y(r)\nu _\beta |\left(𝒱(r_0r)\right)^{}|\nu _e\nu _e|𝒱(r_0)|\nu _\alpha `$. This shows that the perturbary contribution to the transition probability is proportional to $`\mathrm{\Delta }S`$, and therefore, largely independent of the form of the perturbation. Note that perturbation theory only holds for $`\mathrm{\Delta }S`$ small compared to $`rAr`$, i.e., $`\mathrm{\Delta }A/Ar/\lambda `$.

## Acknowledgments

T.O. and W.W. greatly acknowledge the support for the visit at KTH - SCFAB, where a large part of this work was carried out, as well as for the warm hospitality.

This work was supported by the Swedish Foundation for International Cooperation in Research and Higher Education (STINT) \[T.O.\], the Wenner-Gren Foundations \[T.O.\], the “Sonderforschungsbereich 375 für Astro-Teilchenphysik der Deutschen Forschungsgemeinschaft” \[T.O. and W.W.\], and the Swedish Natural Science Research Council (NFR), Contract No. F 650-19981428/2001 \[H.S.\].