Density dependent hadronic models and the relation between neutron stars
  and neutron skin thickness
arXiv:0704.0407v1  [nucl-th]  3 Apr 2007
Density dependent hadronic models and the relation between neutron stars and
neutron skin thickness
S.S. Avancini,1 J.R. Marinelli,1 D.P. Menezes,1 M.M.W. Moraes,1 and C. Providência2
Depto de F́ısica - CFM - Universidade Federal de Santa Catarina Florianópolis - SC - CP. 476 - CEP 88.040 - 900 - Brazil
Centro de F́ısica Teórica - Dep. de F́ısica - Universidade de Coimbra - P-3004 - 516 - Coimbra - Portugal
In the present work we investigate the main differences in the lead neutron skin thickness, binding
energy, surface energy and density profiles obtained with two different density dependent hadron
models. Our results are calculated within the Thomas-Fermi approximation with two different
numerical prescriptions and compared with results obtained with a common parametrization of the
non-linear Walecka model. The neutron skin thickness is a reflex of the equation of state properties.
Hence, a direct correlation between the neutron skin thickness and the slope of the symmetry
energy is found. We show that within the present approximations the asymmetry parameter for low
momentum transfer polarized electron scattering is not sensitive to the model differences.
PACS number(s): 21.65.+f,24.10.Jv,95.30.Tg,26.60.+c
I. INTRODUCTION
The relation between neutron star properties which are
obtained from adequate equations of state (EoS) and the
neutron skin thickness has long been a topic of investi-
gation in the literature. The details of this relation and
the important quantities to be discussed have been well
established in [1], where it was shown that the difference
between the neutron and the proton radii, the neutron
skin thickness, is linearly correlated with the pressure of
neutron matter at sub-nuclear densities. This is so be-
cause the properties of neutron stars are obtained from
appropriate EoS whose symmetry energy depends on the
density and also controls the size of the neutron skin
thickness in heavy and asymmetric nuclei, as 208 Pb, for
instance. It is important to remember that the EoS in
neutron stars is also very isospin asymmetric due to the
β- equilibrium constraint.
Hence, isospin asymmetry plays a major role in the un-
derstanding of the density dependence of the symmetry
energy and the consequences it may arise [2]. In [3, 4]
it was shown that the models that yield smaller neutron
skins in heavy nuclei tend to yield smaller neutron star
radii due to a softer EoS.
Neutron stars are believed to have a solid crust formed
by nonuniform neutron rich matter in β-equilibrium
above a liquid mantle. In the inner crust nuclei coex-
ist with a gas of neutrons which have dripped out. The
properties of this crust as, for instance, its thickness and
pressure at the crust-core interface depend a lot on the
density dependence of the EoS used to describe it [4, 5].
On the other hand, it is well known [6, 7] that the ex-
istence of phase transitions from liquid to gas phases in
asymmetric nuclear matter (ANM) is intrinsically related
with the instability regions which are limited by the spin-
odals. Instabilities in ANM described within relativistic
mean field hadron models, both with constant and den-
sity dependent couplings at zero and finite temperatures
have already been investigated [7] and it was shown that
the main differences occur at finite temperature and large
isospin asymmetry close to the boundary of the insta-
bility regions. In neutral neutron-proton-electron (npe)
matter the electrons are also included. In a thermody-
namical calculation the instabilities almost completely
disappear due to the high electron Fermi energy [8].
However, in a dynamical calculation which includes the
Coulomb interaction and allows for independent neutron,
proton and electron fluctuations [9, 10], it is seen that
the electron dynamics tends to restore the short wave-
length instabilities although moderated by the high elec-
tron Fermi energy.
Moreover, it is also known that the liquid-gas phase
transition in ANM can lead to an isospin distillation phe-
nomenon, characterized by a larger proton fraction in
the liquid phase than in the gas phase. This is due to
the repulsive isovector channel of the nuclear interaction
[11–13].
In a recent work the spinodal section and related quan-
tities, as the neutron to proton density fluctuations re-
sponsible for the distillation effect, has been studied
within different relativistic models [8]. It was shown that
the distillation effect within density dependent relativis-
tic models decreases with density above a nuclear density
of ∼ 0.02−0.03 fm−3, a result similar to the one obtained
with the SLy230a parametrization of Skyrme interaction
[14] and contrary to the results found with the more com-
mon relativistic parametrizations with no density depen-
dent coupling parameters. In the last case the distillation
effect becomes always larger as the density increases.
Also, the behavior of the symmetry energy obtained
with density dependent models is closer to what one ob-
tains with non-relativistic models than with other rel-
ativistic models with constant couplings [7]. In an at-
tempt to understand this behavior, a comparison be-
tween the non-relativistic Skyrme effective force and rel-
ativistic mean field models at subsaturation densities was
performed [15]. It was shown that the relativistic mod-
els could also be reduced to an energy density functional
similar to the one describing the Skyrme interaction.
http://arxiv.org/abs/0704.0407v1
There have already been some efforts in order to com-
pare nuclear matter and finite nuclei properties obtained
both with relativistic and non-relativistic models [16, 17]
but there is no clear or obvious explanations for the differ-
ences. At very low densities both, the relativistic and the
non-relativistic approaches predict a non-homogeneous
phase commonly named pasta phase, formed by a com-
petition between the long-range Coulomb repulsion and
the short-range nuclear attraction [18].
Based on the above arguments, it is very important
that an accurate experimental measurement of the neu-
tron skin thickness is achieved. This depends on a precise
measurement of both the charge and the neutron radius.
The charge radius is already known within a precision of
one percent for most stable nuclei, using the well-known
single-arm and non-polarized elastic electron scattering
technique as well as the spectroscopy of muonic atoms
[19] . For the neutron radius, our present knowledge has
an uncertainty of about 0.2 fm [20]. However, using po-
larized electron beams it is possible to obtain the neutron
distribution in nuclei in a fairly model independent way,
as first discussed in [21] and, as a consequence, to obtain
the desired neutron radius. In fact, the Parity Radius
Experiment (PREX) at the Jefferson Laboratory [22] is
currently running to measure the 208Pb neutron radius
with an accuracy of less than 0.05 fm, using polarized
electron scattering.
In the present work, we use two different hadronic
models that incorporate density dependence in differ-
ent ways. The first one, to which we refer next as the
TW model is a density dependent hadronic model with
the meson-to-nucleon couplings explicitly dependent of
the density [23, 24]. In the following it is used to cal-
culate the neutron skin thickness of 208Pb, which is a
neutron-rich heavy nucleus. This model was chosen be-
cause it is based on a microscopic calculation, fits well
many nuclei properties and, as stated above, has shown
to provide results which are different from the usual NL3
[25] and TM1 [26] parametrizations for the non-linear
Walecka model (NLWM), having a richer density depen-
dence of the symmetry energy than most of the rela-
tivistic nuclear models. The original motivation for the
development of this density dependent hadronic model
[27, 28] was to reproduce results obtained with the rel-
ativistic Dirac-Brueckner Hartree-Fock (DBHF) theory
[29]. Later the DBHF calculations for nuclear matter
were taken only as a guide for a suitable parametrization
of the density dependence of the meson-nucleon coupling
operators [24, 30]. Moreover, density dependent hadronic
models can also be a useful tool in obtaining EoS for
neutron stars even if hyperons are to be considered [32],
which is not the case if NL3 or TM1 are used. Both, NL3
and TM1, can only be used if the EoS is restricted to ac-
commodate neutrons, protons and the leptons necessary
to enforce β-stability. Once hyperons are included, the
nucleons acquire a negative effective mass above∼ 3−4ρ0
densities [33, 34], where ρ0 is the nuclear saturation den-
sity.
The second model, that we refer to as NLωρ model,
includes non-linear σ − ρ and ω − ρ couplings [3, 35–
37] which allow to change the density dependence of the
symmetry energy of the most common parametrizations
of the NLWM that show essentially a linear behavior of
the symmetry energy with density. However, the symme-
try energy determines the behavior of isospin asymmetric
matter and therefore is intrinsically related to the char-
acteristics of the EoS that can describe neutron stars.
Within this model the authors of [3] have shown that
the neutron skin thickness of 208Pb was sensitive to the
isovector channel of the nuclear interaction and there was
a correlation between neutron skin thickness of nuclei and
properties of neutron stars.
For the sake of completeness, the results of the present
work, whenever possible are compared with the results
obtained with the NL3 parametrization of the NLWM,
known to describe finite nuclei properties well.
We perform two different numerical calculations to ob-
tain the 208Pb properties: a Thomas-Fermi approxima-
tion based on the liquid-gas phase transition developed
in [38] and a Thomas-Fermi approximation based on a
method proposed in [39], where a harmonic oscillator ba-
sis is used. We restrict ourselves to the Thomas-Fermi
approximation because, as we show in the Results sec-
tion at the end of the paper, for the purpose of obtaining
correct surface energy and neutron-skin thickness, it is
almost as good as the solution of the Dirac equation.
At this point it is worth mentioning that the scalar-
isovector δ mesons, which play an important role in the
isospin channel, could also be incorporated in our work
as done in [7, 9, 40] but in order to make the compar-
isons among different approximations as simple as possi-
ble, they will be included in a future work. Finally, as
we are interested in nuclei ground state properties, all
calculations are performed at zero temperature.
II. THE TW DENSITY DEPENDENT
HADRONIC MODEL
Next we describe the main quantities of the TWmodel,
which has density dependent coupling parameters. The
Lagrangian density reads:
L = ψ̄
i∂µ − ΓvV µ −
τ · bµ
−e (1 + τi3)
− (M − Γsφ)
(∂µφ∂
µφ−m2sφ2)−
m2vVµV
Bµν ·Bµν +
m2ρbµ ·bµ−
µν (1)
where φ, V µ, bµ and Aµ are the scalar-isoscalar, vector-
isoscalar and vector-isovector meson fields and the pho-
ton field respectively, Ωµν = ∂µVν − ∂νVµ , Bµν =
∂µbν − ∂νbµ − Γρ(bµ × bν), Fµν = ∂µAν − ∂νAµ and
τp3 = 1, and τn3 = −1. The parameters of the model
are: the nucleon mass M = 939 MeV, the masses of the
mesons ms, mv, mρ, the electromagnetic coupling con-
stant e =
4π/137 and the density dependent coupling
constants Γs, Γv and Γρ, which are adjusted in order
to reproduce some of the nuclear matter bulk properties
shown in Table I, using the following parametrization:
Γi(ρ) = Γi(ρsat)hi(x), x = ρ/ρsat, (2)
hi(x) = ai
1 + bi(x+ di)
1 + ci(x + di)2
, i = s, v (3)
hρ(x) = exp[−aρ(x− 1)], (4)
with the values of the parameters mi, Γi(ρsat), ai, bi, ci
and di, i = s, v, ρ given in [24]. This model does not
include self-interaction terms for the meson fields (i.e.
κ = 0, λ = 0 and ξ = 0 ) as in NL3 or TM1 parametriza-
tions for the NLWM.
The field equations of motion follow from the Euler-
Lagrange equations. When they are obtained, some care
has to be taken since the coupling operators depend on
the baryon fields ψ̄ and ψ through the density. When the
partial derivatives of L are performed relatively to the
fields ψ̄ and ψ, they yield extra terms due to the func-
tional dependence of the coupling operators. The new
terms are absent in the usual Quantum Hadrodynamic
(QHD, NLWM) models [25, 26, 31]. The equations of
motion for the fields read:
µ +m2φ)φ = Γsψ̄ψ, (5)
µν +m2vV
µ = Γvψ̄γ
µψ, (6)
µν +m2ρb
ψ̄τγµψ, (7)
ψ̄(1 + τ3)γ
µψ, (8)
[γµ(i∂
µ − Σµ)−M∗]ψ = 0 , (9)
whereM∗ =M−Γsφ. Notice that in the equation of mo-
tion for the baryon field ψ the vector self-energy consists
of two terms, Σµ = Σ
µ + Σ
µ , where:
Σ(0)µ = ΓωVµ +
τ · bµ +
(1 + τ3)Aµ, (10)
ΣRµ =
V νjν +
bν · jν3 −
where Σ
µ is the usual vector self-energy, ρ̂uµ = jµ with
u2 = 1 jν = ψ̄γνψ, j
3 = ψ̄τγ
νψ and, as a result of
the derivative of the Lagrangian with respect to ρ a new
term appears, ΣRµ , which is called rearrangement self-
energy and has been shown to play an essential rôle in
the applications of the theory. This term guarantees the
thermodynamical consistency and the energy-momentum
conservation. For more detailed calculations, at zero and
finite temperatures, please refer to [41].
In the static case there are no currents in the nucleus
and the spatial vector components are zero. Therefore,
the mesonic equations of motion become:
∇2φ = m2sφ− Γsρs, (12)
∇2V0 = m2vV0 − Γvρ, (13)
∇2b0 = m2ρb0 −
ρ3, (14)
∇2A0 = −eρp, (15)
where ρs =< ψ̄ψ > is the scalar density, ρ = ρp + ρn,
ρ3 = ρp − ρn and ρp and ρn are the proton and neutron
densities.
A. Thomas-Fermi approximation
We first define the functional
Ω = E − µpBp − µnBn, (16)
where E is the energy, µp (µn) is the proton (neutron)
chemical potential and Bp (Bn) is the proton (neutron)
number. Within the semi-classical Thomas-Fermi ap-
proximation, the energy of the nuclear system with par-
ticles described by the one-body phase-space distribution
function f(r,p, t) at position r, instant t with momentum
p is given by
(2π)3
fi(r,p, t)
p2 +M∗2 + Vi
(∇φ)2 +m2sφ2 − (∇V0)2 −m2vV 20
−(∇b0)2 −m2ρb20 − (∇A0)2
where
Vp = ΓvV0 +
b0 + eA0 , Vn = ΓvV0 −
γ = 2 refers to the spin multiplicity and the distribu-
tion functions for protons and neutrons are
fi = θ(k
Fi(r) − p2), i = p, n .
In this approach, the scalar, proton and neutron densities
become:
ρs(r) =
i=p,n
∫ kFi(r)
with ǫ =
p2 +M∗2 and
d3rρi, ρi(r) =
k3Fi(r).
From the above expressions we get for (16)
(∇φ)2 − (∇V0)2 − (∇b0)2 − (∇A0)2
+ Vef
Vef =
2 −m2vV 20 −m2ρb20
− µpρp − µnρn
i=p,n
∫ kFi
dpp2ǫ+ΓvV0ρ+Γρ
ρ3 + eA0ρp (18)
Minimization of Ω with respect to kFi(r), i = p, n,
gives rise to the following conditions
k2Fp +M
∗2 − ΓvV0 −
b0 − eA0 − ΣR0
k2Fn +M
∗2 − ΓvV0 +
b0 − ΣR0
where the rearrangement term is
ΣR0 =
ρ V0 +
− ∂ Γs
ρs φ.
From the above equations we obtain kFp = 0 and kFn = 0
or, for kFp or kFn different from zero,
k2Fp +M
∗2 + ΓvV0 +
b0 + eA0 +Σ
0 , (19)
k2Fn +M
∗2 + ΓvV0 −
b0 +Σ
0 . (20)
The values of kFp and kFn are obtained inverting these
two last equations.
Such density dependences in the coupling parameters
do not affect the energy functional but of course affect
its derivative such as the pressure density and the chem-
ical potentials. As already discussed in the literature
[7–9, 32], the rearrangement term is crucial in obtaining
different behaviors in physical properties related to the
chemical potentials or to their derivatives with respect to
the density, such as spinodal regions, as compared with
the more common NL3 or TM1 parametrizations.
III. NLωρ MODEL
The Lagrangian density that incorporates the extra
non-linear σ − ρ and ω − ρ couplings [3, 35–37] reads
L = ψ̄
i∂µ − gvV µ −
τ · bµ
−e (1 + τi3)
− (M − gsφ)
(∂µφ∂
µφ−m2sφ2)−
κφ3 − 1
λφ4 − 1
m2vVµV
µ − 1
Bµν ·Bµν +
m2ρbµ · bµ −
+g2ρbµ · bµ[Λsg2sφ2 + Λvg2vVµV µ], (21)
where Ωµν , Bµν and Fµν are defined after eq.(1). The
parameters of the model are again the masses and the
couplings, which are now constants, i.e., gs replaces Γs,
gv replaces Γv and gρ replaces Γρ. Non-linear σ terms are
also included. We have followed the prescription of [3],
where the starting point was the NL3 parametrization
and the gρ coupling was adjusted for each value of the
coupling Λi studied in such a way that for kF = 1.15
fm−1 (not the saturation point) the symmetry energy is
25.68 MeV. In the present work we set Λs = 0 as in [37].
Notice that other possibilities for this model with σ − ρ
and ω − ρ couplings have already been discussed in the
literature as in [4], for instance.
The mesonic equations of motion in the Thomas-Fermi
approximation become
∇2φ = m2sφ− gsρs +
φ3 (22)
∇2V0 = m2vV0 − gvρ+ 2Λvg2v V0 g2ρb20, (23)
∇2b0 = m2ρb0 −
ρ3 + 2Λvg
0 , (24)
∇2A0 = −eρp, (25)
and the expression for the energy reads
∫ kFi(r)
(2π)3
p2 +M∗2
(∇φ)2 +m2sφ2 − (∇V0)2 −m2vV 20
−(∇b0)2 −m2ρb20 − (∇A0)2
+gvV0ρ+
ρ3b0 + eA0ρp
φ4 − Λvg2vV 20 g2ρb20
. (26)
All other expressions are very similar to the ones ob-
tained from the TW model and can be read off from
them bearing in mind that the density dependent cou-
plings have to be replaced by the constant couplings. In
particular the chemical potentials do not contain the re-
arrangement term ΣR0 .
IV. NUMERICAL RESULT VIA A
NUCLEATION PROCESS
At this point, eqs. (12-15) for the TW model and eqs.
(22-25) for the NLωρmodel have to be solved numerically
in a self-consistent way and hence, initial and boundary
conditions for each equation are necessary. One of the
methods we use here is based on a prescription given in
[38], where these conditions are obtained from a situation
of phase coexistence in a mean field approximation with
classical meson fields and no electromagnetic interaction.
The method is well explained in [38] and, as we are using
different models here, just the main equations are written
next.
For the TW model, the equilibrium equations for ho-
mogeneous matter for the fields are:
m2sφ− Γs ρs = 0, (27)
m2vV0 − Γv ρ = 0, (28)
m2ρb0 −
ρ3 = 0, (29)
and for the energy and pressure density:
E = 1
∫ kFi
p2 +M∗2
V 20 +
b20, (30)
∫ kFi
V 20 +
+ρΣR0 . (31)
For the NLωρ model, the equilibrium equations for ho-
mogenous matter, energy density and pressure become:
m2sφ− gsρs +
φ3 = 0, (32)
m2vV0 − gvρ+ 2Λvg2v V0 g2ρb20 = 0, (33)
m2ρb0 −
ρ3 + 2Λvg
0 = 0, (34)
E = 1
∫ kFi
p2 +M∗2
2 −m2vV 20 −m2ρb20
+ gvV0ρ+
φ4 − Λvg2vV 20 g2ρb20. (35)
∫ kFi
V 20 +
φ3 − λ
φ4 + Λvg
0. (36)
Based on the geometrical construction and Gibbs con-
ditions for phase coexistence, i.e., the pressure and both
chemical potentials are equal in both phases, we build
the binodal section given in Fig. 1. Notice that we have
defined the proton fraction of the system as
. (37)
The binodal section yields the boundary conditions which
we need. For the same pressure, two points, with differ-
ent proton fractions are found. For each of these points,
the meson fields and the densities are well defined and
used as the initial and boundary conditions in eqs. (12-
15), which are then solved. Once the meson fields are
obtained, all the quantities that depend on them, as the
energy, pressure densities, chemical potentials, baryonic
densities, etc are also computed. The solution is a droplet
with a certain proton fraction surrounded by a gas of
neutrons. If stable nuclei are calculated, the gas vanishes
because the energy of the system lies below the neutron
drip line and the finite nuclei properties are easily calcu-
lated. This is the general method, but the results depend
strongly on the model used because of the reasons dis-
cussed in Section VI.
V. NUMERICAL RESULT WITHIN A
HARMONIC OSCILLATOR BASIS
Here a different prescription for solving the equations
of motion and the thermodynamical quantities within the
Thomas-Fermi approximation is used. According to [39],
meson field equations of motion of the Klein-Gordon type
with sources can be carried out by an expansion in a com-
plete set of basis states. The harmonic oscillator func-
tions with orbital angular momentum equal to zero are
then chosen. The oscillator length is given by
, b0 =
, (38)
where M is the nucleon mass and ω0 is the oscillator
frequency. The meson fields and their corresponding in-
homogeneous part can be expanded as
Λ(r) =
ΛnRn0(r), SΛ(r) =
SΛnRn0(r), (39)
where Λ(r) = φ(r), V0(r), b0(r) and
Rnl(r) =
l+1/2
n−1 (x
2)exp(−x2/2), (40)
where x = r/b0 is the radius measured in units of the
oscillator length,
Nnl =
2(n− 1)!/(l + n− 1/2)! (41)
is the normalization constant and Lmn (x
2) are the asso-
ciated Laguerre polynomials. For the calculation of the
meson fields l = 0 in the expressions given below. Once
the ansatz given by eqs.(39) are substituted into eqs.(12-
14), a set of inhomogeneous equations is obtained:
Hnn′Λn′ = SΛn (42)
where
Hnn′ = δnn′
b−2B (2(n− 1) + 3/2) +m
+δnn′+1b
n(n+ 1/2) + δn+1n′b
n′(n′ + 1/2).
vΛ   =0.01
vΛ   =0.025
 0  0.1  0.2  0.3  0.4  0.5
FIG. 1: Binodal section for the NL3, TW and NLωρ
parametrizations.
Only the massive fields can be calculated with this
method because the convergence of the Coulomb field,
which has a long range, is very slow. The Green’s func-
tion method is then chosen to describe the electromag-
netic interaction:
A0(r) = e
r′2dr′ρp(r
′)Gc(r, r
′), (44)
Gc(r, r
1/r for r > r′
1/r′ for r′ > r.
VI. RESULTS
A. Parity Violating Electron Scattering and the
Neutron Radius
We start this section by defining the asymmetry for
polarized electron scattering of a hadronic target as
A = dσ+/dΩ− dσ−/dΩ
dσ+/dΩ+ dσ−/dΩ
, (46)
where dσ±/dΩ is the differential cross section for initially
polarized electrons with positive(+) and negative (−) he-
licities. As the electromagnetic interaction is not sensi-
tive to the above difference, the asymmetry becomes de-
pendent of the weak interaction between the electron and
the target. Moreover, we know from the Standard Model
that the neutral Z-boson couples more strongly to the
neutron than to the proton. Those reasonings were then
used in [21] to first propose a clean way to determine
the neutron distribution in nuclei. If we consider elastic
scattering on an even-even target nucleus, the asymmetry
can be written in the form:
V + β
ρn(q)
ρp(q)
]. (47)
In the above expression, G, α, a and β
V are Standard
Model coupling constants as defined in [21], q is the trans-
ferred momentum by the electron to the nucleus and,
ρn(p)(q) =
d3r j0(qr)ρn(p)(r), (48)
ρn(p)(r) being the neutron (proton) distribution in con-
figuration space and j0 the spherical Bessel function of
order zero. It is then clear that a small q measurement
of the asymmetry gives the neutron radius of the distri-
bution once the proton radius is well known. The proton
and neutron mean-square radius are defined as
R2i =
d3rr2ρi(r)
d3rρi(r)
, i = p, n. (49)
The neutron skin thickness is defined as
θ = Rn −Rp. (50)
In the PREX experiment mentioned in the Introduc-
tion, the asymmetry is expected to be measured at
q ≈ 0.4 fm−1 [22]. Also, because the target is a heavy nu-
cleus (208Pb), the above results for the asymmetry should
be reconsidered for a detailed comparison with the exper-
iment, since they were obtained using a Plane Wave Born
Approximation for the electron [43]. For our present pur-
poses, eq. (47) is sufficient to illustrate the sensitivity to
the different model parametrizations and is used next in
the presentation of our numerical results.
The surface energy per unit area of the droplets in
the small surface thickness approximation, excluding the
electromagnetic field, reads [38]
. (51)
However, as the electromagnetic interaction does not con-
tribute to surface properties directly, we have kept the
same definition for the surface energy.
In Table II we show the neutron and proton radius, the
neutron skin thickness, the binding energy and the sur-
face energy obtained within the Thomas-Fermi approxi-
mation and the two different numerical prescriptions de-
scribed in the previous sections. All the results are sensi-
tive to the numerical calculation although the analytical
approximation is the same. When the nucleation method
is performed, the neutron radius is systematically larger,
what results in a thicker neutron skin. This is correlated
with the fact that the surface energy is lower within the
nucleation calculation than within the harmonic oscilla-
tor method. Within the same numerical prescription,
the neutron skin thickness is smaller with the TW model
than with the NL3. As the coupling strength Λv increases
in the NLωρ model, the results move from the original
NL3 to the TW results for all quantities, except the pro-
ton radius, which oscillates a little. We have also in-
cluded the results obtained with the HS parametrization
[44] because we have used this parametrization in order
to compare the TF and the Dirac results for the cross sec-
tions, as discussed in the following. As this parametriza-
tion is known not to give as good results as the other
parametrizations of the NLWM for finite nuclei, we do
not comment on the results it provides. Notice that the
experimental radius for the protons is obtained from the
charge radius Rc and it is given by Rp =
R2c − 0.64 in
fm [39]. Our results can be compared with experimental
and other theoretical results found in the literature. The
proton radius, which is known to better than 0.001 fm is
better described within the TW model. This quantity is
practically independent of the ω− ρ interaction strength
in the NLωρ model as far as the HO numerical prescrip-
tion is used. The neutron radius, on the other had, is
strongly model dependent with drastic consequences in
the neutron skin thickness calculation. The experimen-
tal values for θ are still very uncertain and all our re-
sults fall inside the experimental confidence interval. We
shall comment on possible restrictions to the neutron skin
thickness in the next section. NL3 provides the best re-
sults for the binding energy.
In [25], the results shown for the proton and neutron
radius are respectively 5.52 and 5.85 fm, yielding a skin
of 0.33 fm, larger than ours. Notice, however, that in
[25] the Dirac equation was explicity solved. In [4], the
authors obtained a value of 0.21 fm for the neutron skin
thickness and a binding energy of -7.89 MeV within a
different parametrization of the NLωρ model. Again in
this case the Dirac equation was solved.
In Fig. 2 we show the difference between neutron and
proton densities at the Pb surface for the models dis-
cussed in the present work with the Thomas-Fermi ap-
proximation solved in a harmonic oscillator basis. While
the curves deviate a little in between 6.0 and 8.0 fm, at
the very surface they are similar, but a small discrepancy,
reflecting the differences in the neutron skin can be seen.
In Fig. 3 we display again the difference between neu-
tron and proton densities within both numerical calcula-
tions of the TW and NL3 models. These two Thomas-
Fermi calculations should have given more similar results.
However the nucleation method predicts a very small sur-
face energy for the NL3 parametrization, and therefore,
a large radius. This may be related to the choice of the
boundary conditions and a deeper comparison between
the two methods will be pursued.
Next we present our results for the asymmetry given
by eq. (47) as a function of the transfered momentum.
We begin with Fig.4 which displays the results for the HS
parametrization of the Walecka model. The curve labeled
no structure means the case where Zρn(r) = Nρp(r) and
the other two curves are obtained within the TF approx-
imation and the full solution of the Dirac equation in
6 7 8 9
0.000
0.005
0.010
0.015
0.020
0.025
 NL ( v )
 NL ( v )
r(fm)
FIG. 2: Difference between neutron and proton densities ob-
tained with the Thomas-Fermi approach solved in a harmonic
oscillator basis for the models discussed in the present work.
the Hartree approximation. At the momentum trans-
fer values of recent experimental interest (around ≃ 0.4
fm−1), the curves are almost identical. A careful analy-
sis of the same results in a different scale shows us that
the asymmetry changes 12 and 11 percent respectively
within the Dirac and TF approximations in comparison
with the no structure case. Since it is the measurement
of the asymmetry in this low momentum transfer region
that will provide the accurate result for the neutron skin
thickness, we have restricted our calculations to the TF
approximation, as stated in the Introduction.
In Fig. 5a we show the asymmetry obtained with the
NL3 model for both numerical calculations in the TF ap-
proximation, i.e, nucleation and HO expansion methods.
In this case, the agreement is very satisfactory even for
larger q-values, although the small numerical discrepan-
cies is reflected in a ∼ 10 percent difference in the pre-
dicted neutron skin thickness, as can be seen from Table
II. Finally, in Fig. 5b our results for the NLωρ (using
two different values for the ω− ρ coupling constant) and
the TWmodels within the HO numerical prescription are
shown. Again, at low momentum transfers, all curves co-
incide. However, it should be noticed that even for two
different model parametrizations which lead us to identi-
cal neutron skin thicknesses, a measurement of the asym-
metry in a higher q-region with a modest experimental
precision, can distinguish between them. Also, we should
expect that the asymmetry presents more structure in
this high momentum transfer region if we solve the Dirac
equation instead of using the TF approach, once the high
q value region is much more sensitive to the central part
of the neutron distribution, which is known to be flat in
the TF approximation. These differences can be seen in
Fig.4.
 NL3 nucl
 TW nucl
 NL3 HO
 TW HO
r(fm)
FIG. 3: Difference between neutron and proton densities ob-
tained with the Thomas-Fermi approach solved with both nu-
merical prescriptions for the TW model.
VII. DIFFERENT EOS, DIFFERENT NEUTRON
SKINS
For the sake of completeness, at this point, we discuss
some of the differences between the TW, the NLωρ mod-
els and the NL3 parametrization of the NLWM.
From Fig. 1 one can see that the largest possible pres-
sure for a phase coexistence in the TW model is much
lower, and appears at a lower proton fraction than the
NL3 model. This gives rise to a thinner crust within
the TW model, which may imply that the more exotic
pasta shapes will not form [5]. The NLωρ model goes
on a different direction, i.e., the pressure becomes higher
than the one obtained with the NL3 as the Λv coupling
is turned on.
Although the nuclear matter properties fitted to
parametrize the models are quite similar (see Table I),
the way the EoS behaves when extrapolated to higher
or lower densities can vary a lot from a density depen-
dent hadron model to one of the parametrizations of the
NLWM. Moreover, as seen from Table I, although the ef-
fective mass at saturation density is lower with the TW
than with the NL3, it can accommodate hyperons if an
EoS for stellar matter is necessary, contrary to the usual
0.0 0.5 1.0 1.5 2.0
 no structure
 HS-Dirac 
 HS-TF A
q(fm-1)
0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50
2.5x10-7
5.0x10-7
7.5x10-7
1.0x10-6
1.3x10-6
1.5x10-6
 no structure
 HS-Dirac 
 HS-TF 
q(fm-1)
FIG. 4: Parametrization HS, comparison Thomas-Fermi-HO
versus Dirac-HO
NL3 parametrization [32–34].
TABLE I: Nuclear matter properties.
NL3 NLωρ TW
[25] [36] [24]
Λv = 0.01 Λv = 0.02 Λv = 0.025
B/A (MeV) 16.3 16.3 16.3 16.3 16.3
ρ0 (fm
−3) 0.148 0.148 0.148 0.148 0.153
K (MeV) 271 271 271 271 240
Esym. (MeV) 37.4 34.9 33.1 32.3 32.0
M∗/M 0.60 0.60 0.60 0.60 0.56
L (MeV) 118 88 68 61 55
Ksym (MeV) 100 -46 -53 -34 -124
Another quantity of interest in asymmetric nuclear
matter is the nuclear bulk symmetry energy, shown in
Table I for the saturation point. The differences in the
symmetry energy at densities larger than the nuclear sat-
uration density is still not well established, but has al-
ready been extensively discussed in the literature even
0.0 0.5 1.0 1.5 2.0
q(fm-1)
 NL3 nucl
 NL3 HO
 no structure
0.0 0.5 1.0 1.5 2.0
q(fm-1)
 NL ( v=0.01)
 NL ( v=0.025)
FIG. 5: Asymmetry obtained with a) NL3 with both numer-
ical prescriptions and b)parametrizations NLωρ and TW
for the TW model [7, 8, 16, 32]. Again, for the sake
of completeness we reproduce these results here because
the neutron skin thickness and the neutron star EoS are
related by this quantity [1–4], which is usually defined
as Esym = 12
∂2E/ρ
, with δ = −ρ3/ρ = 1 − 2yp. The
symmetry energy can be analytically rewritten as
Esym =
ρ, (52)
for the TW model and as
Esym =
ρ (53)
with the effective ρ-meson mass defined as [3]
= m2ρ + 2g
for the NLωρ model. In both cases
kFp = kF (1 + δ)
1/3, kFn = kF (1− δ)1/3,
with kF = (1.5π
2ρ)1/3 and ǫF =
k2F +M
∗2. In equa-
tions (52) and (53) the second term dominates at large
densities. It is seen that the non-linear ρ − ω terms
introduce a non-linear density behavior in the symme-
try energy of the NLWM parametrizations such as NL3
and TM1. In TW the non-linear density behavior en-
ters through the density dependent coupling parameters.
These non-linear density behavior is important because
the linear behavior of NL3 and TM1 parametrizations
predicts too high symmetry energy at densities of impor-
tance for neutron star matter which has direct influence
on the proton fraction dependence with density. From
Fig. 6, it is easily seen that the symmetry energy ob-
tained with the TW model behaves in a very different
way, as compared with NL3. In [4] a relation between
the symmetry energy and the nuclear binding energy is
discussed : the harder the EoS, the more the symmetry
(fm  )−3ρ
vΛ   =0.01
vΛ   =0.025
 0  0.05  0.1  0.15  0.2  0.25  0.3  0.35  0.4
FIG. 6: Symmetry energy for the NL3, TW and NLωρmodels.
energy rises with density. The density dependence dis-
cussed in [4] is of the type introduced in [3, 36] through
the inclusion of a σ− ρ and/or ω− ρ couplings and then,
similar with the NLωρ model discussed here. One can
observe that as the strength of the coupling increases,
the symmetry energy gets closer to the TW curve. In
fact, in [8] it was shown that once this kind of coupling
is introduced with a reasonable strength, the symmetry
energy at low densities tends to behave as the TW model.
The symmetry energy can be expanded around the nu-
clear saturation density and reads
Esym(ρ) = Esym(ρ0) +
ρ− ρ0
ρ− ρ0
where L and Ksym are respectively the slope and the
curvature of the nuclear symmetry energy at ρ0 and they
are calculated from
L = 3ρ0
∂Esym(ρ)
|ρ=ρ0 Ksym = 9ρ20
∂2Esym(ρ)
|ρ=ρ0 .
These two quantities can provide important information
on the symmetry energy at both high and low densi-
ties because they characterize the density dependence of
the energy symmetry. In a recent work [49], the authors
found a correlation between the slope of the symmetry
energy and the neutron skin thickness. In their work 21
sets of the non-relativistic Skyrme potential were inves-
tigated and only 4 of them were shown to have L values
consistent with the values extracted from experimental
isospin diffusion data from heavy ion collisions. In fact,
the extracted value was L = 88 ± 25 MeV [50], which
gives a very strong constraint on the density dependence
of the nuclear symmetry energy and consequently on the
EoS as well. A detailed analysis of Table I shows that, if
this constraint is to be taken seriously, neither the NL3
nor the TW model satisfy it. Nevertheless, the NLωρ
slope interpolates beautifully between the NL3 and TW
slope values. Once again it is seen that the increase in
Λv approximates the NL3 model values for the slope and
energy symmetry to the TW values. Moreover, we have
also tried to find a correlation between the θ values shown
in Table II and L values displayed in Table I. We found
that, as far as some numerical imprecision are consid-
ered, larger values of L correspond to larger values of the
neutron skin, as seen in Fig. 7.
Let’s now go back to the problem of solving the dif-
ferential equations within the nucleation numerical pre-
scription. As we need boundary conditions arising from
the liquid-gas phase coexistence in order to solve eqs.
(12-15) for the TW model and eqs. (22-25) for the NLωρ
model, the binodal sections are essential and the spinodal
sections, which separate the regions of stable to unsta-
ble matter are also of interest. If we had displayed the
binodals in a ρp versus ρn plot, as it is done with the
spinodals in Fig 8, we could see that the spinodals sur-
faces lie inside the binodal sections and share the critical
point corresponding to the highest pressure.
In Fig. 8 the spinodals for the three different mod-
els discussed in this work are shown. Once again, some
of these results can also be found in the recent litera-
ture [7, 8], but we include them here to make a direct
link with the binodals. The instability of the ANM sys-
tem is essentially determined by density fluctuations in
the isoscalar channel. Although the spinodals are, by
themselves, not relevant in calculations performed at the
thermodynamical equilibrium, the isospin channel is very
sensitive to the instabilities occurring below the nuclear
saturation density. The spinodal is determined by the
values of pressure, proton fraction and density for which
the determinant of
Fij =
∂ρi∂ρj
, (56)
where F is the free energy density, goes to zero. A de-
tailed analysis of this quantity can be found in [8, 42].
From Fig. 8, it is seen that the instability region in the
ρp/ρn plane, defined by the inner section of the spinodal
curve is larger for the TW than for the NL3 model. The
size of the instability region depends on the derivative
of the chemical potentials with respect to the neutron
and proton densities. At low densities different models
exhibit different behaviors.
The presence of the rearrangement term in the TW
model also plays a decisive role. Even though a rela-
tively large compensation exists between scalar and vec-
tor mesons in the isoscalar channels within the rearrange-
ment term at low densities, the spinodal region is defined
by the derivative of the chemical potential and therefore
of the rearrangement term.
Next we examine the spinodals obtained with differ-
ent coupling strengths for the NLωρ model. As seen in
Fig. 8, there is almost no difference between the different
curves. They all fall around the original NL3 curve but
once again, they tend to the TW curve as the coupling
strength increases. However, contrary to the TW model,
it was shown in [9] that the direction of the instability in
Λ  =0.01
Λ  =0.02
Λ  =0.025v 0.16
 0.18
 0.22
 0.24
 50  60  70  80  90  100  110  120
L (MeV)
FIG. 7: Correlation between the neutron skin θ and the slope
of the symmetry energy L.
NLωρ increases distillation as the density increases, and
the larger the coupling Λv the larger the effect.
Finally, to end this section, let’s make our points clear:
we have used a simple mean field theory approach to ob-
tain the boundary conditions for the equations of motion
of the meson fields in the nucleation prescription. These
boundary conditions depend on the model used and are
intrinsically related with the liquid-gas phase transition
which, in turn, can be well understood by studying the
coexistence surfaces of the corresponding models. On
the other hand, the neutron skin thickness shows a lin-
ear correlation with the slope of the symmetry energy,
as already pointed out in [49] for non-relativistic mod-
els. Based on the different behaviors found with density
dependent hadronic models and the NLWM, an obvious
consequence is the fact that the neutron skin thickness
depends on the choice of the model.
VIII. CONCLUSIONS
We have calculated the 208Pb neutron skin thickness
with two different density dependent hadronic models,
the TW and the NLωρ model, and one of the most used
parametrizations of the NLWM, the NL3. The calcu-
lations were done within the Thomas-Fermi approxima-
tion, which gives quite accurate results for the asymme-
try in the momentum transfer range of interest for the
calculation of neutron skins. In implementing the nu-
merical results two different prescriptions were used: the
first one based on the nucleation process and the second
one based on the harmonic oscillator basis method. We
have seen that when the nucleation method is performed,
the neutron radius is systematically larger, what results
vΛ   =0.01
vΛ   =0.025
ρ (fm   )−3
 0.02
 0.04
 0.06
 0.08
 0  0.02  0.04  0.06  0.08
FIG. 8: Spinodal section in terms of ρp versus ρn for the NL3,
TW and NLωρ models.
in a thicker neutron skin. This is a consequence of the
fact that the surface energy is lower within the nucleation
calculation than within the harmonic oscillator method.
Within the same numerical prescription, the neutron skin
thickness is smaller with the TW model than with the
NL3. As the coupling strength Λv increases in the NLωρ
model, the neutron skin thickness moves from the orig-
inal NL3 towards the TW results. We have also found
that although the neutron skin thickness is model depen-
dent, the asymmetry at low momentum transfers (below
0.5 fm−1) is very similar for all models and all numerical
prescriptions. As q increases, the asymmetry also be-
comes model dependent. The density profiles obtained
from the solution of the Dirac equation exhibits oscil-
lations near the center of the nucleus, behavior which is
not reproduced within the Thomas-Fermi approximation.
This fact shows up in the asymmetry at large momentum
transfers and therefore all the calculations should be re-
produced by solving the Dirac equation. This calculation
is already under investigation.
It is worth mentioning that the neutron skin thickness
has shown to give hints on the equations of state that
are suitable to describe neutron stars. Moreover, in [49]
a correlation between the slope of the symmetry energy
and the neutron skin thickness was found for Skyrme-
type models. We have observed that this correlation was
also present in the density dependent models we have
studied in the present work.
ACKNOWLEDGMENTS
This work was partially supported by CNPq(Brazil),
CAPES(Brazil)/GRICES (Portugal) under project
100/03 and FEDER/FCT (Portugal) under the projects
POCTI/FP/63419/2005 and POCTI/FP/63918/2005.
[1] S. Typel and B.A. Brown, Phys. Rev. C 64, 027302
(2001).
[2] A.W. Steiner, M. Prakash, J.M. Lattimer and P.J. Ellis,
Phys. Rep. 411, 325 (2005).
[3] C.J. Horowitz and J.Piekarewicz, Phys. Rev. Lett. 86,
5647 (2001).
[4] J.Piekarewicz, nucl-th/0607039. Proceedings of the ”In-
ternational Conference on Current Problems in Nuclear
Physics and Atomic Energy” (May 29 - June 3, 2006)
Kyiv, UKRAINE.
[5] F. Duchoin and Haensel, Phys. Lett. B 485, 107 (2000).
[6] Ph. Chomaz, C. Colonna and J. Randrup, Phys. Rep.
389, 263 (2004).
[7] S.S. Avancini, L. Brito, D. P. Menezes and C.
Providência, Phys. Rev. C 70, 015203 (2004).
[8] S.S. Avancini, L. Brito, Ph. Chomaz, D. P. Menezes and
C. Providência, Phys. Rev. C 74, 024317 (2006).
[9] C. Providência, L. Brito, S.S. Avancini, D. P. Menezes
and Ph. Chomaz, Phys. Rev. C 73, 025805 (2006).
[10] L. Brito, C. Providência, A.M.S. Santos, S.S. Avancini,
D. P. Menezes and Ph. Chomaz. Phys. Rev. C (2006),
C 74, 045801 (2006); C. Providência, L. Brito, A.M.S.
Santos, D.P. Menezes and S.S. Avancini, Phys. Rev. C
74, 045802 (2006).
[11] Ph. Chomaz and F. Gulminelli, Phys. Lett. B447, 221
(1999) 221; H. S. Xu, et al, Phys. Rev. Lett. 85, 716
(2000).
[12] C. Ducoin, Ph. Chomaz and F. Gulminelli, Nucl. Phys.
771, 68 (2006).
[13] Ph. Chomaz, Nucl. Phys. A 685, 274c (2001).
[14] E. Chabanat, P. Bonche, P. Haensel, J. Meyer and R.
Schaeffer, Nucl. Phys. A 627, 710 (1997).
[15] C. Providência, D. P. Menezes, L. Brito and Ph. Chomaz,
in preparation.
[16] B.A. Li, C.M. Ko and W. Bauer, Inter. J. Mod. Phys. E
7, 147 (1998).
[17] K. Pomorski, P. Ring, G.A. Lalazissis, A. Baran, Z. Lo-
jewski, B. Nerlo-Pomorska, M. Warda, Nucl. Phys. A
624, 349 (1997).
[18] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys.
Rev. Lett. 50, 2066 (1983); M. Hashimoto, H. Seki, and
M. Yamada, Prog. Theor. Phys.71, 320 (1984).
[19] H. de Vries, C.W. de Jager and C. de Vries, Atomic and
Nuclear Data Tables 36, 495 (1987).
[20] C.J. Horowitz, S.J. Pollock, P.A. Souder and R. Michaels,
Phys. Rev. C 63, 025501 (2001).
[21] T.W. Donnelly, J. Dubach and I. Sick, Nucl. Phys. A503
589 (1989).
[22] K.A. Aniol et al. (HAPPEX) (2005), nucl-
ex/0506010; ibidem, nucl-ex/0506011; R. Michaels,
P.A. Souder and G.M. Urciuoli (2005), URL
http://hallaweb.jlab.org/parity/prex.
[23] H. Lenske and C. Fuchs, Phys. Lett. B 345, 355 (1995);
C. Fuchs, H. Lenske and H.H. Wolter, Phys. Rev. C 52,
3043 (1995).
[24] S. Typel and H. H. Wolter, Nucl. Phys. A656, 331
(1999).
[25] G. A. Lalazissis, J. König and P. Ring, Phys. Rev. C 55,
540 (1997).
[26] K. Sumiyoshi, H. Kuwabara, H. Toki, Nucl. Phys. A 581,
725 (1995).
[27] C. Fuchs, H. Lenske and H.H. Wolter, Phys. Rev. C 52,
3043 (1995).
[28] H. Lenske and C. Fuchs, Phys. Lett. B 345, 355 (1995).
[29] B. ter Haar and R. Malfliet, Phys. Rep. 149, 207 (1987).
[30] T. Niks̆ić, D. Vretenar, P. Finelli and P. Ring, Phys. Rev.
C 66, 024303 (2002).
[31] B. Serot and J.D. Walecka, Advances in Nuclear Physics
16, Plenum-Press, (1986) 1.
[32] S.S. Avancini and D.P. Menezes, Phys. Rev. C 74,
015201 (2006).
[33] D.P. Menezes and C. Providência, Phys. Rev. C 68,
035804 (2003); Braz. J. Phys. 34, 724 (2004).
[34] A.M.S. Santos and D.P. Menezes, Phys. Rev. C 69,
045803 (2004).
[35] C.J. Horowitz and J.Piekarewicz, Phys. Rev.C 64,
062802 (2001).
[36] J.K. Bunta and S. Gmuca, Phys. Rev. C 68, 054318
(2003).
[37] J.K. Bunta and S. Gmuca, Phys. Rev. C 70, 054309
(2004).
[38] D.P. Menezes and C. Providência, Nucl. Phys. A 650,
283 (1999); D.P. Menezes and C. Providência, Phys. Rev.
C 60, 024313 (1999); D.P. Menezes and C. Providência,
Phys. Rev. C 64, 044306 (2001).
[39] Y.K. Gambhir, P. Ring and A. Thimet, Ann. Phys. 198,
132 (1990).
[40] T. Gaitanos, M. Di Toro, S. Typel, V. Baran, C. Fuchs,
V. Greco and H. H. Wolter, Nucl. Phys. A 732, 24
(2004).
[41] S.S. Avancini, M.E. Bracco, M. Chiapparini and D.P.
Menezes, J. Phys. G 30, 27 (2004); S.S. Avancini, M.E.
Bracco, M. Chiapparini and D.P. Menezes, Phys. Rev. C
67, 024301 (2003).
[42] J. Margueron and P. Chomaz, Phys. Rev. C 67, 041602
(2003).
[43] C.J. Horowitz, Phys. Rev. C57, 3430 (1998).
[44] C.J. Horowitz and B.D. Serot, Nucl. Phys. A 368, 503
(1981).
[45] G.Fricke, C. Bernhardt, K.Heilig, L.A. Schaller, L.
Schellinberg, E.B. Shera,C.W. de Jager, At. Data Nucl.
Data Tables 60 (1995)177.
[46] G. Audi, A.H. Waptra, C. Thibault, Nucl. Phys. A 729,
337 (2003).
[47] A. Krasznahorkay et a., Nucl. Phys. A 731, 224 (2004).
[48] V.E. Starodubsky, N.M. Hintz,Phys. Rev.
C49,2118(1994).
[49] L. Chen, C.M. Ko and B. Li, nucl-th/0610057.
[50] M.B. Tsang et al., Phys. Rev. Lett. 92, 062701 (2004).
TABLE II: 208 Pb properties
model approximation Rn Rp θ B/A σ
(fm) (fm) (fm) MeV Mev/fm2
NL3 TF+nucleation 5.88 5.65 0.24 -7.77 0.76
NL3 TF+HO 5.79 5.57 0.22 -7.79 0.96
NLωρ, Λv = 0.01 TF+HO 5.77 5.57 0.20 -7.73 0.98
NLωρ, Λv = 0.02 TF+HO 5.75 5.57 0.17 -7.65 0.99
NLωρ, Λv = 0.025 TF+HO 5.74 5.58 0.16 -7.63 1.00
TW TF+nucleation 5.71 5.50 0.22 -6.42 1.08
TW TF+HO 5.68 5.52 0.16 -7.46 1.10
HS TF+HO 5.70 5.47 0.24 -6.10 1.37
exp.[45] 5.44
exp. [46] -7.87
exp. [47] 0.12± 0.07
exp. [48] 0.20± 0.04
