ar
X
iv
:n
uc
l-t
h/
04
01
02
4v
1
1
2
Ja
n
20
04
Mean field and pairing properties in the crust of neutron stars
F. Montani
Institut d’ Astronomie et d’ Astrophysique,
Universite Libre de Bruxelles,
Campus de la Plaine CP-226,
Boulevard du Triomphe, B-1000 Brussels, Belgium
C. May and H. Mu¨ther
Institut fu¨r Theoretische Physik,
Universita¨t Tu¨bingen, D-72076 Tu¨bingen, Germany
Properties of the matter in the inner crust of a neutron star are investigated in a Hartree-Fock
plus BCS approximation employing schematic effective forces of the type of the Skyrme forces.
Special attention is paid to differences between a homogenous and inhomogeneous description of
the matter distribution. For that purpose self-consistent Hartree Fock calculations are performed
in a spherical Wigner-Seitz cell. The results are compared to predictions of corresponding Thomas
Fermi calculations. The influence of the shell structure on the formation of pairing correlations in
inhomogeneous matter are discussed.
I. INTRODUCTION
The determination of the equation of state (EoS) for nucleon matter is an important ingredient for various investiga-
tions of astrophysical objects. A lot of attention has been paid to the EoS of baryonic matter at densities above the
saturation density of nuclear matter (ρ0 ≈ 0.16 fm−3). These densities should be relevant to describe the interior of
neutron stars and are of particular interest as they give rise for speculations about exotic phases of baryonic matter
including the existence of kaon condensates or quark matter[1, 2, 3].
However, it is not only this regime of very high densities also the crust of neutron stars should be a very intriguing
phase of baryonic matter. At those lower densities and low temperatures the free energy should be reduced by allowing
for a phase of inhomogeneous baryonic matter. This nonuniform matter should consist out of a lattice of quasi-nuclei
embedded in a gas of electrons and possibly also a sea of neutrons. Thomas Fermi calculations, which are based on
non-relativistic or relativistic mean field calculations for the homogenous system, exhibit these features[4, 5, 6, 7, 8].
The presence of neutron superfluidity in the crust of neutron stars seems to be well established. This is based on
investigations solving the BCS gap equation in homogenous nuclear matter[9, 10, 11, 12]. The precise knowledge of
the pairing gap is very important for the determination of the cooling rate of neutron stars. The superfluidity of
the material of the crust should have a significant influence on the rotation frequencies of the star. The formation of
glitches may be related to the vortex pinning of the superfluid phase in the inhomogeneous crust.
The aim of the investigations presented in this paper is to explore basic properties of this inhomogeneous phase of
baryonic matter on a level beyond the Thomas Fermi approximation. What are the consequences of shell effects in
the quasi-nuclei in the systems for the bulk-properties of this material? How do they influence the proton fraction in
the β equilibrium? What are the consequences of these quasi-nuclei on the pairing properties of the system?
To answer this kind of questions we perform nuclear structure calculations in a Wigner-Seitz cell of spherical shape.
The assumption of such a spherical cell is not designed to allow for inhomogeneities of the matter which is different
from the formation of quasi-nuclei or bubbles in the homogenous matter. Therefore we do not consider a formation
of rods or slabs of increased density[8]. Also at this stage we do not aim at a complete survey of this matter covering
a large region of various densities. Instead we wish to explore some features of nuclear structure calculations beyond
the Thomas Fermi approach for a few selected examples.
For that purpose we perform Hartree-Fock and Hartree-Fock plus BCS calculations in an appropriate basis of the
spherical Wigner-Seitz cell. For the NN interaction we use rather simple effective forces of the Skyrme type[13, 14].
Special attention is paid to the pairing gaps resulting from calculations of the homogenous matter and the non-
homogeneous system.
After this introduction we discuss in section 2 the techniques and results of Hartree-Fock calculations for the β-stable
matter within a spherical Wigner-Seitz cell. The determination of pairing properties is presented in section 3, which
is followed by a short summary and conclusions in section 4.

2II. HARTREE-FOCK CALCULATIONS IN A SPHERICAL BOX
The single-particle wave functions for the nucleon considered in our calculations are expanded in a complete basis of
orthonormal states defined within a spherical box of radius R. Such a orthonormal set of basis functions, which are
regular at the origin in the center of the box is given by
Φiljm(~r) = 〈~r|iljm〉 = Ril(r)Yljm(ϑ, ϕ) . (1)
In this equation Yljm represent the spherical harmonics including the spin degrees of freedom by coupling the orbital
angular momentum l with the spin to a single-particle angular momentum j. The radial wave functions Ril are given
by the spherical Bessel functions, Ril(r) ∼ jl(kir) for the discrete momenta ki, which fulfill
Ril(R) = Niljl(kiR) = 0 . (2)
The normalization constant
Nil =
{ √2√
R3jl−1(kiR)
for l > 0 ,
ipi
√
2√
R3
for l = 0 ,
(3)
ensures that the basis functions are orthogonal and normalized within the box,
∫ R
0
d3rΦ∗iljm(~r)Φi′l′j′m′(~r) = δii′δll′δjj′δmm′ . (4)
Adopting this basis of eigenstates for the kinetic energy one can try to describe homogenous nuclear matter with a
constant density ρpi and ρν for protons and neutrons, respectively, by filling all basis states with momenta ki below
the corresponding Fermi momentum kFpi and kFν . The local density is then given by
ρ(r) =
∑
ilj
Θ(kF − ki)(2j + 1)R2il(r) , (5)
where Θ stands for the Heaviside function. Results for this local density considering a spherical box with radius R=25
fm and trying to describe a system with a density of 0.06 fm−3 is displayed in Fig. 1. From this figure one can see that
the local density is in reasonable agreement with the mean value except at the border of the cavity at r = R. Since
all wave functions by construction are bound to disappear at this borderline, the same is also true for the density.
In order to cure this deficiency one could try an alternative basis for the radial functions by the boundary condition
that the radial derivative vanishes at the surface of the box
R˜il(r) = N˜iljl(k˜ir) with
∂R˜il
∂r (R) = 0 . (6)
So all these basis function exhibit an extremum at the surface of the Wigner-Seitz cell which leads to a maximum of
the local density at r = R if these functions are employed in (5) as on can see from the curve labeled “Basis 2” in
Fig. 1.
Bonche and Vautherin[15] suggested to use a mixed basis by employing e.g basis states with the boundary condition
(2) for states with even l and those with the boundary condition (6) for states with odd orbital angular momentum l.
This recipe cures the deviation of the local density from the mean value at the surface (see also Fig. 1). Nevertheless
the representation of the homogenous phase of nuclear matter within a spherical box of finite size leads to slight
fluctuations in the local density. Furthermore one must be aware of shell effects, which are due to the finite size of
the spherical box considered.
In order to explore the influence of these shell effects on the evaluation of expectation values, we have considered
spherical boxes of various radii and calculated the binding energy per nucleon for homogenous nuclear matter using
the representation of the plane wave states discussed above. We have considered values for the radius R ranging from
15 fm to 25 fm. Most of the results turned out to be rather insensitive to the choice of the basis ((2), (6) or mixture).
Therefore, if not stated differently, we will present results for the basis (2) only.
As a first example we consider the homogenous system of neutrons, protons and electrons in β equilibrium. This
means that for any value of the baryon density considered
ρ = ρpi + ρν , (7)

3we determine the proton abundance
Ypi =
ρpi
ρ , (8)
by the requirement that the Fermi energy for the neutrons is identical to the Fermi energy of the protons plus the
Fermi energy of the electrons with the density of electrons being identical to the density of protons.
Results for the energy per nucleon and the proton abundance are displayed in figures 2 and 3, respectively. For
the range of densities considered the results obtained in the spherical boxes of different size agree rather well with
the corresponding observables calculated in the infinite system. A significant discrepancy is only observed in the
calculation of the total energy at higher densities considering boxes with small radii. In these cases the surface effects
displayed in Fig.1 as well as the shell effects lead to energies which are to small as compared to the result for the
infinite system.
After this test of the box basis to describe the homogenous system, we now turn to the Hartree-Fock (HF) description
within such a spherical box. The HF single-particle wave functions are expanded in the basis of (1) by
Ψαljm(~r) =
N∑
i=1
cαiljΦiljm(~r) . (9)
The number of basis states N is chosen to guarantee that the results are not affected by this limitation. Results of
self-consistent HF calculations are displayed in Fig. 4 showing the density of protons and neutrons as a function of
the radial distance from the center of a box with radius R=21 fm. the particle numbers for protons (Z=20 in this
example) and neutrons (N=388) have been determined to fulfill the condition of β stability, employing the Fermi
energies evaluated in the HF calculation.
The density distribution for the protons is different from zero only near the center of the Wigner Seitz cell. In this
region also the density of the neutrons is considerably larger than in the rest of the spherical box. Therefore we
can interpret this configuration as a quasi-nucleus embedded in a neutron see. The single-particle energy spectrum
for the protons exhibit a clear shell structure up to energies around the Fermi energy. The density of states around
the corresponding Fermi energy is larger for the neutrons. In fact it is close to the one of the homogeneous matter
described in the spherical box of this size.
For the example of Z=20 protons and N=388 neutrons in a Wigner Seitz cell with radius R=21 fm considered in
Fig. 4 we obtain a global density of ρ = 0.0105 fm−3 and a abundance for the protons of Yp = 0.049. The energy per
nucleon for this self-consistent solution of the HF equations is about 2 MeV per nucleon below the energy which was
obtained for the homogenous distribution of matter in β equilibrium.
Fig. 5 exhibits this gain in energy due to the formation of quasi-nuclei in β-stable nuclear matter. For small densities
around ρ = 0.01 fm−3 the HF solution with localized quasi-nuclei yields an energy which is about 3 MeV per nucleon
below the energy of the corresponding homogenous matter. This gain in energy decreases with increasing density to
around 1 MeV per nucleon at ρ = 0.04 fm−3. The results are rather independent on the size of the Wigner Seitz cell
under consideration.
The Hartree-Fock calculations allowing for an inhomogeneous distribution of matter, however, also yield different
results for the proton abundances as compared to the results obtained for homogenous nuclear matter in β equilibrium
(see Fig. 6). Thomas Fermi calculation, which allow for an inhomogeneous distribution of matter, show proton
abundances which are slightly above those derived from the calculation of homogenous matter. This indicates that
the fluctuations in the density, which are taken into account in the Thomas Fermi calculations, tend to enhance the
proton abundances, as regions with high densities contain a larger fraction of protons. However, the main effect in the
enhancement of the proton abundances observed in the Hartree Fock approach is not described by the Thomas Fermi
model. It originates from the pronounced shell structure of the proton single-particle energies. The Hartree-Fock
calculation favors distributions of matter with a quasi-nucleus showing a closed shell for the protons.
III. PAIRING IN THE CRUST OF NEUTRON STARS
After we have noticed the effects of shell structure on the decomposition of the nuclear material, we would like to
address in this section the question how this shell structure can affect the pairing properties of the nuclear material.
In particular we consider neutron-neutron pairing for neutron pairs with total momentum equal to zero in 1S0 partial
wave for the relative motion. Using the standard BCS approach for homogenous matter the pairing gap ∆(k) for a

4pair of neutrons with relative momentum k = |~k| is obtained by solving the gap equation[9, 11, 12]
∆(k) = − 2π
∫ ∞
0
dk′ k′2 V (k, k′) ∆(k
′)
2
√
(ε′k − εF )
2 +∆(k′)2
. (10)
Here V (k, k′) denotes the matrix elements of the NN interaction in the 1S0 partial wave, εk the single-particle energy
for a nucleon with momentum k and εF the Fermi energy. Instead of using the matrix elements for a realistic NN
interaction which is fitted to the NN scattering data, we have chosen to employ an effective interaction which is of
zero range, density dependent and has been proposed by Bertsch and Esbensen [16, 17]
V (r1, r2) = v0
[
1− η
(
ρ((r1)
ρ0
)α]
δ(r1 − r2) , (11)
where v0, η and α are adjustable parameters and ρ0 is the saturation density of nuclear matter. For such a zero-
range interaction a cut-off parameter εC must be introduced in the gap equation to truncate the integral to states
with energy εk less than εC . These four parameters can be constrained to reproduce the NN scattering length and
Garrido et al.[18] determined various sets of parameters which reproduce pairing gaps for nuclear and neutron matter
calculated with realistic NN interactions. We used the parameters v0 = 481 MeV fm3, η = 0.45, α = 0.47 and εC =
60 MeV and verified that this set of parameters leeds to pairing gaps for homogenous neutron matter, which are in
fair agreement with those evaluated for the CD Bonn interaction[19].
We now turn to the solution of the HF+BCS equations[20], which in the spherical Wigner-Seitz cell have the form
(εnlj − εF )unlj +∆nljvnlj = Enljunlj
− (εnlj − ǫF ) vnlj +∆nljunlj = Enljvnlj , (12)
where
Enlj =
√
(εnlj − εF )2 +∆2nlj ,
is the energy of the quasiparticle state with quantum numbers n, l, j and εnlj the corresponding single-particle energy
determined in the HF equations. From the coefficients unlj , vnlj and the corresponding single-particle wave functions
Ψnljm (see eq.(9)) one can calculate the anomalous density
χ(r) =
∑
nlj
(2j + 1)unljvnlj
2
Ψ2nlj(r) , (13)
and the expression for the normal density is changed into
ρ(r) =
∑
nlj
(2j + 1)v2nljΨ2nlj(r) . (14)
Using an interaction of zero range like (11) leads to a definition of a local gap-function
∆(r) = V (r)χ(r) , (15)
from which one can calculate the state-dependent pairing gaps
∆nlj =
∫
∆(r)Ψ2nlj(r) r2 dr . (16)
The BCS equations (12) - (16) have to be solved together with the HF equations in a self-consistent way. The Fermi
energy εF is adjusted to reproduce the requested particle numbers or densities for protons and neutrons. These
equations can be reduced to an approach, which we will call BCS with plane waves in the spherical box, by restricting
the radial wave functions for the single-particle states to the spherical Bessel functions defined in (2) or (6).
Results for the gap function ∆(r) defined in (15) for β-stable matter of a density ρ = 0.02 fm−3 are displayed in Fig. 7.
The gap function evaluated for the plane wave basis are presented by the dashed line. We find that this gap-function
fluctuates around 1.8 MeV with a sharp drop at the boundary of the spherical box. This drop is of course related to
the boundary condition (2) of all wave functions. As already discussed for the local density above, this deficiency can

5be cured by using a mixed basis, switching between (2) and (6) for states with even and odd parity (see lower panel
of Fig. 7).
Solving the BCS plus Hartree-Fock equations one obtains a gap-function with values which are suppressed by about
25 percent in the area of the quasi-nuclei. From this result one may speculate that the inhomogeneous matter leads to
regions of enhanced densities in which the formation of pairing correlations is suppressed to some extent. Therefore
these regions might be considered as nuclei for the formation of normal matter within the superfluid phase of neutron
matter. Also one could imagine that these quasi-nuclei could lead to vortex pinning in the rotation of the superfluid
crust of neutron stars.
Before drawing any conclusions of this kind, we would like to compare the results for the state-dependent gaps ∆nlj
displayed in Figs. 8 and 9 for the plane wave plus BCS and HF plus BCS calculations, respectively. Comparing these
results, one finds that plane wave plus BCS calculation leads to a gap ∆nlj which is almost independent of the state.
The calculation of the inhomogeneous system on the other hand leads to a few single-particle states, which are much
more bound than the lowest states in the homogeneous calculation (-42 MeV as compared to -9.5 MeV). These deeply
bound states are localized near the origin. This large binding of the single-particle states which leads to a reduction
of the density of states is also responsible for the reduction of the pairing gap for these low-lying states. On the other
hand, the states with single-particle energies close to the Fermi energy εF show a value for the gap parameter which
is very close to the one derived for the homogenous system.
This implies that the reduction of the gap-function ∆(r) in the region of the quasi-nuclei discussed above manifests
itself mainly in a reduction of the pairing gap for the deeply bound single-particle states. The evaluation of response
functions and other observables, however, is in general sensitive to the gap at the Fermi energy. Therefore one cannot
expect major differences in the behavior of the inhomogeneous as compared to the homogenous system of matter with
respect to the pairing properties of the neutrons. This shall be different for smaller densities, where the Fermi energy
for the neutrons drops to values at which the shell structure of the single-particle energies is significant.
IV. CONCLUSIONS
The transition from homogenous to inhomogeneous matter as it occurs in the crust of neutron stars has been in-
vestigated. Special attention has been paid to the consequences of shell effects, which occurring assuming a phase
of quasi-nuclei embedded in a sea of neutrons. For that purpose Hartree-Fock plus BCS calculations are performed
assuming a basis of single-particle states for a spherical Wigner-Seitz cell. It is observed that the shell effects lead to
a substantial increase for the proton abundance in β-stable baryonic matter.
Shell effects are also responsible for a decrease of the localized pairing gap in the region of the quasi-nuclei. A more
detailed analysis, however, shows that this reduction of the local gap function ∆(r) is mainly due to a reduction of
the pairing gap for the deeply bound single-particle states. The pairing properties for the single-particle states close
to the Fermi energy are similar to those obtained for a homogenous description of the system.
The present studies are based on simple parameterizations for the effective NN interaction. It may be of interest to
study whether the predictions for the transition from homogenous to inhomogeneous matter also hold, when more
realistic NN interactions are employed.
The use of a spherical Wigner-Seitz cell is a source of various deficiencies. The boundary conditions at the border of
the cell lead to fluctuations in the density, which complicates the comparison with the infinite homogenous system.
This could be improved by employing a Cartesian basis, which is more involved from the numerical point of view.
The use of a Cartesian basis, however, would also allow the study of non-spherical structures.
We would like to thank Stephane Goriely, Mike Pearson and Mathieu Samyn for very helpful discussions. One of us
(F.M.) has been supported by a stipend of the “Deutsche Akademische Austauschdienst” (DAAD, Germany). Also
this financial support is gratefully acknowledged
[1] M. Prakash, I. Bombaci, M. Prakash, P. J. Ellis, J. M. Lattimer, and R. Knorren, Phys. Rep. 280, 1 (1997).
[2] H. Heiselberg and M. Hjorth-Jensen, Phys. Rep. 328, 237 (2000).
[3] S. Pal, M. Hanauske, I. Zakout, H. Stocker, and W. Greiner, Phys. Rev. C 60, 015802 (1999).
[4] K. Oyamatsu, Nucl. Phys. A561, 431 (1993).
[5] K. Sumiyoshi, K. Oyamatsu, and H. Toki, Nucl. Phys. A595, 327 (1995).
[6] H. Shen, Phys. Rev. C 65, 035802 (2002).
[7] M. Pi, X. Vinas, M. Barranco, A. Perez-Canyellas, and A. Polls, Astronomy Astrophysics Suppl 64, 439.
[8] G. Watanabe, K. Iida, and K. Sato, Nucl. Phys. A676, 455 (2000).

6[9] M. Baldo, J. Cugnon, A. Lejeune, and U. Lombardo, Nucl. Phys. A515, 409 (1990).
[10] V.A. Kodel, V.V. Kodel, and J.W. Clark, Nucl. Phys. A598, 390 (1996).
[11] Ø. Elgarøy and M. Hjorth-Jensen, Phys. Rev. C57, 1174 (1998).
[12] J. Kuckei, F. Montani, H. Mu¨ther, and A. Sedrakian, Nucl. Phys. A723, 32 (2003).
[13] T.H.R. Skyrme, Nucl. Phys. 9, 615 (1959).
[14] D. Vautherin and D.M. Brink, Phys. Rev. C5, 626 (1972).
[15] P. Bonche and D. Vautherin, Nucl. Phys. A 372, 496 (1981).
[16] G.F. Bertsch and H. Esbensen, Ann. Phys. (N.Y.) 209, 327 (1991).
[17] H. Esbensen, G.F. Bertsch and K. Hencken, Phys.Rev. C 56, 3054 (1997).
[18] E. Garrido, P. Sarriguren, E. Moya de Guerra, and P. Schuck, Phys. Rev. C 60 064312 (1999).
[19] R. Machleidt, Phys. Rev. C63, 024001 (2001).
[20] P. Ring and P. Schuck, The Nuclear Many Body Problem (Springer Verlag, New York, 1980).

70 5 10 15 20 25
r [fm]
0
0.05
0.1
D
en
sit
y
ρ
[fm
-
3 ]
Basis 1
Basis mixed
Basis 2
FIG. 1: The local density for homogeneous matter (ρ = 0.06 fm−3) calculated according to (5) using Basis 1 (2), Basis 2 (6)
and the mixed basis as described in the text
0 0.02 0.04 0.06 0.08
density ρ [fm-3]
3
4
5
6
7
8
9
10
En
er
gy
/N
uc
le
on
[M
eV
]
R=15
R=17
R=19
R=21
infinite matter
FIG. 2: Energy per nucleon for homogenous matter in β equilibrium. Results for infinite matter are compared to those obtained
in spherical boxes of different radii. The Skyrme I has been used for the NN interaction.

80 0.02 0.04 0.06 0.08
density ρ [fm-3]
0
0.01
0.02
0.03
0.04
0.05
0.06
pr
ot
on
a
bu
nd
an
ce
Y
p
infinite matter
R=15
R=17
R=19
R=21
FIG. 3: Proton abundance Yp for homogenous matter in β equilibrium. Further comments see Fig. 2.
0 5 10 15 20
r [fm]
0
0.02
0.04
0.06
0.08
0.1
de
ns
ity
ρ

[fm
-
3 ]
Protons, Z=20
Neutrons, N=388
FIG. 4: Density distributions derived from Hartree Fock calculations for nuclear matter in β equilibrium.

90 0.01 0.02 0.03 0.04 0.05 0.06
density ρ [fm-3]
0
2
4
6
8
En
er
gy
/N
uc
le
on
[M
eV
]
R=15
R=17
R=19
R=21
inhomogeneous
homogeneous
FIG. 5: Energy per nucleon of β-stable nuclear matter in the homogenous phase (symbols connected by dashed line) and
inhomogeneous phase (symbols connected by solid line).
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
density ρ [fm-3]
0
0.05
Pr
ot
on
a
bu
nd
an
ce
Y
p
R=15
R=17
R=19
R=21
infinite matter
Thomas Fermi
FIG. 6: Proton abundances of β-stable nuclear matter in the homogenous case of homogenous infinite matter (solid line),
inhomogeneous matter determined in the Thomas Fermi approximation and in Hartree-Fock calculations employing Wigner
Seitz cells of different radii.