# Nuclear incompressibility and sound speed in uniform matter and finite
nuclei

G. Grams<EMAIL_ADDRESS>Univ Lyon, Univ Claude Bernard Lyon 1,
CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France  Institut
d’Astronomie et d’Astrophysique, CP-226, Université Libre de Bruxelles, 1050
Brussels, Belgium  R. Somasundaram<EMAIL_ADDRESS>Univ Lyon,
Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622,
Villeurbanne, France  J. Margueron Univ Lyon, Univ Claude Bernard Lyon 1,
CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France  E. Khan
IJCLab, Université Paris-Saclay, CNRS/IN2P3, 91405 Orsay Cedex, France

###### Abstract

We have extended the compressible liquid-drop model (CLDM) with a density-
dependent surface term (eCLDM), which allows for a unified description of both
the nuclear ground state energies and the incompressibility modulus in finite
nuclei $K_{A}$. We analyse the role of the nuclear empirical parameters, e.g.,
$K_{\mathrm{sat}}$, $Q_{\mathrm{sat}}$, $L_{\mathrm{sym}}$ and
$K_{\mathrm{sym}}$, which contribute to the bulk properties, as well as the
role of the finite size contributions. For the bulk properties, the density
and isospin dependencies of the nuclear incompressibility in infinite matter
are characterized by introducing new empirical parameters, and two new
constraints for the value of $K_{\mathrm{sym}}$ are suggested. For finite
nuclei, we employ a Bayesian approach coupled to a Markov-Chain Monte-Carlo
(MCMC) exploration of the parameter space to confront the model predictions of
$K_{A}$ in Zr, Sn and Pb isotopes to the experimental data. We show that
$Q_{\mathrm{sat}}\approx-950\pm 200$ MeV describes the experimental
measurements of $K_{A}$ in these isotopes. This value is different from the
ones deduced from phenomenological nuclear energy density functionals,
suggesting a possible explanation of their difficulty to accurately describe
Zr, Sn and Pb data all together. In addition we explore the impact of a
fictitious measurement of the Giant Monopole Resonance energy in 132Sn. We
show that this measurement, provided it is accurate enough, will allow to
better determine $K_{\mathrm{sym}}$ and $K_{\tau}$. Finally we explore the
properties of the sound speed around saturation density and show the important
role of finite size terms in finite nuclei since they reduce the sound speed
to approximately half compared to nuclear matter.

##  I Introduction

The response of nuclear matter to compression and expansion plays a very
important role in many phenomena in nature, from finite nuclei van der Woude
and Harakeh (2001), which can be viewed as non-uniform pieces of nuclear
matter squeezed by the effects of the surface terms, up to astrophysical
nuclear systems such as neutron stars, supernovae or kilonovae Rezzolla  _et
al._ (2018), where nuclear matter explores densities and isospin asymmetries
in extreme regimes. In finite nuclei, the repulsive surface tension and the
Coulomb interaction counter balance the attractive bulk nuclear force and
allow the exploration of densities close to the saturation density of nuclear
matter ($n_{\mathrm{sat}}\approx 0.155$ fm-3 Margueron  _et al._ (2018a)),
while in compact stars, the bulk nuclear force resists gravity for densities
corresponding to several times saturation density. In these examples, the
equilibrium states of these systems represent a balance between the bulk
properties and the action of external forces (finite size terms or
gravitational force). It is then important to quantify precisely the response
of bulk nuclear matter (incompressibility) from analyses of finite nuclei
properties (giant monopole resonances), which is the scope of the present
study.

At first order, the energy required to compress matter from its equilibrium
state, is given by the incompressibility modulus $K_{\mathrm{sat}}$ with an
isospin asymmetry dependence driven by the parameter $K_{\tau}$ Blaizot
(1980): $K_{\mathrm{sat}}+\delta^{2}K_{\tau}$. These nuclear empirical
parameters could be extracted from the analysis of the isoscalar giant
monopole resonance (ISGMR), excited by the scattering of alpha particles, see
for instance Ref. Garg and Coló (2018) and references therein. The relation
between the energy of the ISGMR, $E_{\mathrm{ISGMR}}$, and the
incompressibility in finite nuclei $K_{A}$ is Garg and Coló (2018),

$E_{\mathrm{ISGMR}}=\hbar\sqrt{\frac{K_{A}}{m_{N}\langle
r^{2}\rangle}}\simeq\hbar\sqrt{\frac{5K_{A}}{3m_{N}R_{A}}}\,,$ (1)

where $m_{N}$ is the nucleon mass and $\langle r^{2}\rangle$ is the mean
square radius of the density distribution in finite nuclei. The last
expression is obtained assuming a flat density distribution up to $R_{A}$, as
in the compressible liquid-drop model (CLDM). Considering a leptodermous
expansion as in the liquid-drop model, the incompressibility modulus in finite
nuclei $K_{A}$ can be expressed as Blaizot (1980)

$K_{A}=K_{\mathrm{sat}}+K_{\tau}\delta^{2}+K_{\mathrm{Coul}}\frac{Z^{2}}{A^{4/3}}+K_{\mathrm{surf}}A^{-1/3}+\dots$
(2)

where $K_{\mathrm{sat}}$ and $K_{\tau}$ are the bulk contributions which we
aim to extract from experimental data, $K_{\mathrm{Coul}}$ is the Coulomb
repulsive contribution and $K_{\mathrm{surf}}$ the surface attractive
contribution. This leptodermous expansion is however difficult to employ for
the determination of $K_{\mathrm{sat}}$ and $K_{\tau}$ from experimental
measurements of $E_{\mathrm{ISGMR}}$, as Eqs. (1) and (2) may suggest,
especially since the term $K_{\mathrm{surf}}$ is difficult to fix from the few
existing experimental data. The situation is different in the case of the
leptodermous expansion of the nuclear mass, since more than 2000 nuclei have
been measured Huang  _et al._ (2021). In a recent analysis Li and Xie (2021),
precise values for $K_{\mathrm{sat}}$ and $K_{\tau}$ have been obtained from
Eq. (2), fixing $K_{\mathrm{surf}}=cK_{\mathrm{sat}}$ (with $c\approx-1.2\pm
0.12$ Patra  _et al._ (2002)) and $K_{\mathrm{Coul}}\approx-5.2\pm 0.7$ MeV
Sagawa  _et al._ (2007). This indicates that the key quantity which would
allow the use of such an empirical relation is the surface term. In this
paper, we investigate the impact of the surface term in the CLDM framework and
its role to reproduce experimental data.

By using the energy density functional approach, the first precise extraction
of $K_{\mathrm{sat}}$ gave $K_{\mathrm{sat}}=210\pm 30$ MeV Blaizot (1980),
corrected to $240\pm 20$ MeV later on as a good compromise between data in
208Pb and 90Zr Garg and Coló (2018). The isospin dependence of the
incompressibility $K_{\tau}$ is more difficult to determine from experimental
data. It has recently been extracted from systematical exploration of Sn
isotopic chain, giving $K_{\tau}\approx-550\pm 100$ MeV Garg and Coló (2018).
However, there are several unsettled questions: by using non-magic nuclei, the
analysis of the data requires the understanding of many-body correlations
(pairing, deformation, etc.) on the incompressibility of finite nuclei. The
question of the isoscalar and isovector properties of the incompressibility is
also important since the density and the isospin asymmetry distributions in
finite nuclei are different from one to another. In addition, a systematic
difference between the incompressibility extracted from 208Pb and from 120Sn
(which tends toward $K_{\mathrm{sat}}\simeq 205$ MeVAvogadro and Bertulani
(2013)) remains which origin is still not well understood. This issue could be
of similar origin as the systematical dispersion for $K_{\mathrm{sat}}$
obtained by using different models, where $K_{\mathrm{sat}}\approx 220$ MeV is
preferred by Gogny forces Goriely  _et al._ (2009) while non-linear
relativistic mean field models favor $K_{\mathrm{sat}}\approx 250$ MeV
Lalazissis  _et al._ (2005). It was suggested that these systematical
differences could be related to the different density dependence of the
models, encoded in the nuclear empirical parameter (NEP) $Q_{\mathrm{sat}}$
Khan  _et al._ (2012); Khan and Margueron (2013). It is indeed a general
result that a large uncertainty on a high order NEP impacts the precise
determination of lower order ones Margueron and Gulminelli (2019).

In uniform matter (UM), the incompressibility $K_{\mathrm{UM}}$ is defined as
the second derivative of the energy density
$\epsilon_{\mathrm{UM}}=E_{\mathrm{UM}}/V$ as

$\displaystyle K_{\mathrm{UM}}(n,\delta)$ $\displaystyle=$ $\displaystyle
9n\frac{\partial^{2}\epsilon_{\mathrm{UM}}(n,\delta)}{\partial n^{2}}\,,$ (3)
$\displaystyle=$
$\displaystyle\frac{18}{n}P_{\mathrm{UM}}(n,\delta)+9n^{2}\frac{\partial^{2}e_{\mathrm{UM}}(n,\delta)}{\partial
n^{2}}\,,$ (4)

where $n$ is the isoscalar density $n=n_{n}+n_{p}$ and $\delta$ the isovector
parameter $\delta=(n_{n}-n_{p})/n$, the energy per particle is
$e_{\mathrm{UM}}=\epsilon_{\mathrm{UM}}/n$, and the pressure $P_{\mathrm{UM}}$
is defined as

$P_{\mathrm{UM}}=n^{2}\frac{\partial e_{\mathrm{UM}}}{\partial n}\,.$ (5)

Note that $K_{\mathrm{UM}}=K_{\mathrm{sat}}$ if $n=n_{\mathrm{sat}}$ and
$\delta=0$. In finite nuclei, the isovector parameter is noted
$\delta_{A}=(N-Z)/A$.

In the absence of external forces, such as gravity for instance, matter
minimizes its energy (mechanical equilibrium) by imposing $P_{\mathrm{UM}}=0$.
We note $n_{\mathrm{eq}}^{\mathrm{UM}}$ the equilibrium density in symmetric
(SM) and isospin asymmetric (AM) matter. The latter always deals with small
isospin asymmetries $|\delta_{A}|\lesssim 0.3$ as expected in finite nuclei.
In nuclear matter and at equilibrium, the first term in Eq. (4) vanishes but
in finite nuclei however the equilibrium density $n_{\mathrm{eq}}^{A}$ is
slightly different from the one in uniform matter
$n_{\mathrm{eq}}^{\mathrm{UM}}$, due to the presence of finite size terms
which contribute to the pressure. This effect shifts
$n_{\mathrm{eq}}^{\mathrm{UM}}$ by about 10% at maximum and impacts the value
of the bulk incompressibility in finite nuclei. One could then view the finite
size terms as an ”external” force probing the response of the bulk.
Consequently, there is a contribution of the finite size terms to the
incompressibility in finite nuclei, in addition to the density and isospin
asymmetry dependence of the bulk term Blaizot (1980). In addition, the
equilibrium density $n_{\mathrm{eq}}^{A}$ in finite nuclei varies around
$n_{\mathrm{eq}}^{\mathrm{UM}}$ through the nuclear chart, modifying the value
of the energy in the bulk. Since this value is controlled at first order by
the incompressibility modulus, the energy of finite nuclei in their ground
state also contains a contribution originating from the incompressibility of
nuclear matter, in addition to the symmetry energy and to the finite size
terms. This contribution is difficult to extract from microscopic approaches,
e.g., energy density functional, shell model approaches, as well as ab-initio
ones, but it could be more visible in macroscopic models such the CLDM that we
employ in this study. The fact that the fluctuations in $n_{\mathrm{eq}}^{A}$
impact both the energy $e_{A}$ and the incompressibility $K_{A}$ requires to
employ a model which could describe these two quantities in a unique
framework. This is the motivation for the development of the eCLDM that we
present in this paper.

The CLDM has been shown relevant to describe nuclear masses Myers and
Swiatecki (1969); Weiss and Cameron (1969) and was also employed to study the
clusterized matter present on neutron star crusts Baym  _et al._ (1971);
Steiner (2008); Carreau  _et al._ (2019); Grams  _et al._ (2022a, b). Many
variations of the model can be found in the literature, however, it has been
argued by Blaizot Blaizot (1980) that the CLDM is not appropriate to
accurately extract the incompressibility modulus $K_{\mathrm{sat}}$ from
finite nuclei. The reason lies in the contribution of the density dependent
surface term to the incompressibility, which is absent in most of the
macroscopic models. In the present work, we however construct an extended CLDM
(eCLDM) with a density-dependent surface tension allowing to describe both
nuclear masses and incompressibilities. Furthermore, the bulk term of the
present model is described with the meta-model Margueron  _et al._ (2018a), an
energy density functional in which the parameters of the model are the
empirical parameters of nuclear matter. The meta-model has the advantage of
being flexible enough to allow an independent variation of the NEP and can
thus be used to easily perform a sensitivity analysis of the individual impact
of the NEP on the incompressibility $K_{A}$, as well as extensive searches of
the best parameter sets reproducing experimental data.

The paper is organized as follow: In Sec. II we explore the incompressibility
modulus in nuclear matter in terms of the NEP, or equivalently as a function
of the density and the isospin asymmetry. A new constraint on
$K_{\mathrm{sym}}$ is derived and compared to other existing ones. Following
the line suggested by Blaizot Blaizot (1980), we then address finite nuclei in
Sec. III as described by our eCLDM model (with a density dependent surface
tension), which allows to reproduce both finite nuclei and incompressibility
modulus from the same approach. In Sect. IV we compare the predictions of the
eCLDM to experimental data and analyse the role of the NEP $K_{\mathrm{sat}}$,
$Q_{\mathrm{sat}}$, $L_{\mathrm{sym}}$ and $K_{\mathrm{sym}}$ in a Bayesian
framework. Finally, in Sec. V we discuss the sound speed in both uniform
matter and finite nuclei.

##  II Uniform matter

In this section, we briefly summarize the present understanding of uniform
matter and show how the knowledge of the NEP could be used to explore its
properties around saturation density. We also present an alternative
representation where the reference density is taken to be
$n_{\mathrm{eq}}^{UM}$, the equilibrium density which is a function of
$\delta$, instead of the saturation density $n_{\mathrm{sat}}$ in the usual
approach.

###  II.1 Representation of the nuclear matter properties in terms of the
nuclear empirical parameters

The NEP, e.g., $E_{\mathrm{sat}}$, $E_{\mathrm{sym}}$, are defined as the
coefficients of the series expansion of the energy per particle in SM
($e_{\mathrm{SM}}$) and of the symmetry energy ($e_{\mathrm{sym}}$) as,

$\displaystyle e_{\mathrm{SM}}(n)$ $\displaystyle=$ $\displaystyle
E_{\mathrm{sat}}+\frac{1}{2}K_{\mathrm{sat}}x^{2}+\frac{1}{6}Q_{\mathrm{sat}}x^{3}$
(6) $\displaystyle\hskip
28.45274pt+\frac{1}{24}Z_{\mathrm{sat}}x^{4}+\dots\,,$ $\displaystyle
e_{\mathrm{sym}}(n)$ $\displaystyle=$ $\displaystyle
E_{\mathrm{sym}}+L_{\mathrm{sym}}x+\frac{1}{2}K_{\mathrm{sym}}x^{2}+\frac{1}{6}Q_{\mathrm{sym}}x^{3}$
(7) $\displaystyle\hskip
28.45274pt+\frac{1}{24}Z_{\mathrm{sym}}x^{4}+\dots\,,$

where $x=(n-n_{\mathrm{sat}})/3n_{\mathrm{sat}}$, with $n_{\mathrm{sat}}$
being the saturation density of nuclear matter ($n_{\mathrm{sat}}=0.155\pm
0.005$ fm-3, see for instance Ref. Margueron  _et al._ (2018a)). Note that
choosing $n_{\mathrm{sat}}$ as the reference density for the parameter $x$, is
arbitrary: in Sec. II.2 for instance, we explore another reference density. It
should also be noted that in Eq. (7), the symmetry energy is defined as the
difference between neutron matter (NM) and SM energies, as
$e_{\mathrm{sym}}(n)=e_{\mathrm{NM}}(n)-e_{\mathrm{SM}}(n)$. It can be
expanded in terms of $\delta^{2}$ as
$e_{\mathrm{sym}}(n)=e_{\mathrm{sym},2}(n)\delta^{2}+e_{\mathrm{NQ}}$, where
$e_{\mathrm{sym},2}$ and $e_{\mathrm{NQ}}$ subsume the quadratic and non-
quadratic (NQ) contributions respectively.

It was suggested in Ref. Margueron  _et al._ (2018a) to consider the series
expansion up to order 4 in the density parameter $x$ in order to represent
accurately the energy per particle, the pressure and the sound speed of
existing models up to about $4n_{\mathrm{sat}}$. We adopt this prescription
here as well, even if we do not explore such high densities.

Note that since asymmetric matter is mostly quadratic in $\delta$, as it is
expected to be Somasundaram  _et al._ (2021a), Eqs. (6)-(7) could also be
written in a more compact way,

$\displaystyle e_{\mathrm{UM}}(x,\delta)$ $\displaystyle\approx$
$\displaystyle e_{\mathrm{SM}}(n)+e_{\mathrm{sym}}(n)\delta^{2}\,,$ (9)
$\displaystyle\approx$ $\displaystyle
E(\delta)+L_{\mathrm{sym}}x\delta^{2}+\frac{1}{2}K(\delta)x^{2}+\frac{1}{6}Q(\delta)x^{3}$
$\displaystyle\hskip 56.9055pt+\frac{1}{24}Z(\delta)x^{4}+\dots\,,$

where

$\displaystyle E(\delta)$ $\displaystyle\equiv$ $\displaystyle
E_{\mathrm{sat}}+E_{\mathrm{sym}}\delta^{2}\,,\hskip 8.5359ptK(\delta)\equiv
K_{\mathrm{sat}}+K_{\mathrm{sym}}\delta^{2}\,,$ (10) $\displaystyle Q(\delta)$
$\displaystyle\equiv$ $\displaystyle
Q_{\mathrm{sat}}+Q_{\mathrm{sym}}\delta^{2}\,,\hskip 8.5359ptZ(\delta)\equiv
Z_{\mathrm{sat}}+Z_{\mathrm{sym}}\delta^{2}\,.$ (11)

It should be noted that the above expression of K($\delta)$ is by no means the
true isospin dependence of the incompressibility, as it will be discussed
below. In particular, it neglects the contribution of the pressure which is
different from zero as one gets farther from saturation. It solely represents
the second order term in the density expansion of the energy per particle.

From Eq. (9), one could deduce a similar expression for the energy density
$\epsilon_{\mathrm{UM}}=(1+3x)e_{\mathrm{UM}}n_{\mathrm{sat}}$ as,

$\displaystyle\epsilon_{\mathrm{UM}}(x,\delta)/n_{\mathrm{sat}}$
$\displaystyle=$ $\displaystyle
E(\delta)+L^{\epsilon}(\delta)x+\frac{1}{2}K^{\epsilon}(\delta)x^{2}+\dots\,,$
(12)

where

$\displaystyle L^{\epsilon}(\delta)$ $\displaystyle\equiv$ $\displaystyle
3E_{\mathrm{sat}}+(3E_{\mathrm{sym}}+L_{\mathrm{sym}})\delta^{2}\,,$ (13)
$\displaystyle K^{\epsilon}(\delta)$ $\displaystyle\equiv$ $\displaystyle
K_{\mathrm{sat}}+K_{\mathrm{sym}}^{\epsilon}\delta^{2}\,,$ (14)

with

$\displaystyle K^{\epsilon}_{\mathrm{sym}}$ $\displaystyle\equiv$
$\displaystyle K_{\mathrm{sym}}+6L_{\mathrm{sym}}\,.$ (15)

The $\delta$-dependence of the energy density curvature $K^{\epsilon}(\delta)$
is different from that of the energy per particle curvature $K(\delta)$.
Consequences will be discussed in the following, especially for the
incompressibility modulus in asymmetric matter. It will be shown that
$K^{\epsilon}(\delta)$ do correspond to the isospin dependence of the
incompressibility around saturation density, contrarily to Ksym, which is only
a parameter useful in the expansion (9)

Model | BSK14 | BSK16 | F0 | LNS5 | RATP | SGII | SKI2 | SKO | SLy5
---|---|---|---|---|---|---|---|---|---
Ref. | Goriely  _et al._ (2007) | Chamel  _et al._ (2008) | Lesinski  _et al._ (2006) | Cao  _et al._ (2006) | Rayet  _et al._ (1982) | Van Giai and Sagawa (1981) | Reinhard and Flocard (1995) | Reinhard  _et al._ (1999) | Chabanat  _et al._ (1998)
$E_{\mathrm{sat}}$ (MeV) | -15.85 | -16.05 | -16.03 | -15.56 | -16.05 | -15.59 | -15.76 | -15.83 | -15.98
$n_{\mathrm{sat}}$ (fm-3) | 0.159 | 0.159 | 0.162 | 0.160 | 0.160 | 0.158 | 0.158 | 0.161 | 0.160
$K_{\mathrm{sat}}$ (MeV) | 239 | 242 | 230 | 240 | 240 | 215 | 241 | 223 | 230
$Q_{\mathrm{sat}}$ (MeV) | -359 | -364 | -405 | -316 | -350 | -381 | -339 | -393 | -364
$E_{\mathrm{sym}}$ (MeV) | 30.00 | 30.00 | 32.00 | 29.15 | 29.26 | 26.83 | 33.37 | 31.97 | 32.03
$L_{\mathrm{sym}}$ (MeV) | 43.9 | 34.9 | 42.4 | 50.9 | 32.4 | 37.6 | 104.3 | 79.1 | 48.3
$K_{\mathrm{sym}}$ (MeV) | -152 | -187 | -113 | -119 | -191 | -146 | 71 | -43 | -112
$Q_{\mathrm{sym}}$ (MeV) | 389 | 462 | 658 | 286 | 440 | 330 | 52 | 131 | 501
Table 1: Nuclear empirical parameters for the Skyrme interactions used in the
present work.

The general expressions for the pressure (5) and the incompressibility modulus
(4) in AM could be expressed in terms of the parameter $x$ as,

$\displaystyle P_{\mathrm{UM}}(x,\delta)$ $\displaystyle=$
$\displaystyle\frac{n_{\mathrm{sat}}}{3}(1+3x)^{2}\frac{\partial
e_{\mathrm{UM}}(x,\delta)}{\partial x}\,,$ (16) $\displaystyle
K_{\mathrm{UM}}(x,\delta)$ $\displaystyle=$ $\displaystyle
6(1+3x)\frac{\partial e_{\mathrm{UM}}(x,\delta)}{\partial
x}+(1+3x)^{2}\frac{\partial^{2}e_{\mathrm{UM}}(x,\delta)}{\partial x^{2}}\,.$

Injecting Eq. (9) into the expression for the pressure (16), we obtain

$\displaystyle P_{\mathrm{UM}}(x,\delta)$ $\displaystyle=$
$\displaystyle\frac{n_{\mathrm{sat}}}{3}\left[L_{\mathrm{sym}}\delta^{2}+K^{p}(\delta)x+\frac{1}{2}Q^{p}(\delta)x^{2}\right]$
(18) $\displaystyle\hskip 56.9055pt+o(x^{3})\,,$

where $K^{p}=K^{\epsilon}$ and $Q^{p}$ reads

$\displaystyle Q^{p}(\delta)\equiv
Q^{p}_{\mathrm{sat}}+Q^{p}_{\mathrm{sym}}\delta^{2}\,,$ (19)

with

$\displaystyle Q^{p}_{\mathrm{sat}}$ $\displaystyle\equiv$ $\displaystyle
Q_{\mathrm{sat}}+12K_{\mathrm{sat}}\,$ (20) $\displaystyle
Q^{p}_{\mathrm{sym}}$ $\displaystyle\equiv$ $\displaystyle
Q_{\mathrm{sym}}+18L_{\mathrm{sym}}+12K_{\mathrm{sym}}\,.$ (21)

Note that in finite nuclei, $|\delta_{A}|<0.3$ and densities are explored from
about $2/3n_{\mathrm{sat}}$ up to $n_{\mathrm{sat}}$, which implies
$|x_{A}|\lesssim 0.1$. In finite nuclei, we could therefore perform an
expansion at the same level in $\delta^{2}$ and in $x$.

The equilibrium density in AM is given by the density for which the mechanical
stability is satisfied: $\partial e_{\mathrm{UM}}(x,\delta)/\partial x=0$.
From the expression of the pressure (16) truncated at order $x$, one can
deduce in AM Piekarewicz and Centelles (2009),

$x_{\mathrm{eq}}^{\mathrm{UM}}\approx-\frac{L_{\mathrm{sym}}}{K(\delta)}\delta^{2}\approx-\frac{L_{\mathrm{sym}}}{K_{\mathrm{sat}}}\delta^{2}\,.$
(22)

The equilibrium density is a function of the isospin asymmetry parameter
$\delta$, and it satisfies the limit $n_{\mathrm{eq}}^{\mathrm{UM}}\rightarrow
n_{\mathrm{sat}}$ for $\delta\rightarrow 0$. At order $\delta^{2}$ and $x$,
one obtains for the equilibrium density $n_{\mathrm{eq}}^{\mathrm{UM}}$ in
asymmetric matter,
$n_{\mathrm{eq}}^{\mathrm{UM}}=n_{\mathrm{sat}}[1-3(L_{\mathrm{sym}}/K_{\mathrm{sat}})\delta^{2}]$.

In finite nuclei, the situation is more complex than previously described
since: i) the equilibrium density is different from $n_{\mathrm{sat}}$, due to
the finite size terms and ii) the density is not uniform allowing for surface
contributions to be sizeable. In uniform matter however, only isospin
asymmetry contributes to the shift of the equilibrium density from
$n_{\mathrm{sat}}$, as shown in Eq. (22). While neglecting the contribution of
the finite-size (FS) terms, expression (22) provides a good estimation of the
average densities in finite nuclei Papakonstantinou  _et al._ (2013). In the
next section, this density is named $n_{\mathrm{cl}}$ in the CLDM and we have
$n_{\mathrm{cl}}\approx n_{\mathrm{eq}}^{\mathrm{UM}}$ for large $A$.

The pressure could be decomposed into a SM and an isospin asymmetry terms:

$P_{\mathrm{UM}}=P_{\mathrm{SM}}+P_{\mathrm{sym}}\delta^{2}\,,$ (23)

with

$\displaystyle P_{\mathrm{SM}}(n)$ $\displaystyle=$
$\displaystyle\frac{n_{\mathrm{sat}}}{3}\Big{[}K_{\mathrm{sat}}x+\frac{1}{2}Q^{p}_{\mathrm{sat}}x^{2}+...\Big{]}\,,$
(24) $\displaystyle P_{\mathrm{sym}}(n)$ $\displaystyle=$
$\displaystyle\frac{n_{\mathrm{sat}}}{3}\Big{[}L_{\mathrm{sym}}+K^{p}_{\mathrm{sym}}x+\frac{1}{2}Q^{p}_{\mathrm{sym}}x^{2}+...\Big{]}\,.$
(25)

We have for instance
$P_{\mathrm{sym}}(n_{\mathrm{sat}})=n_{\mathrm{sat}}L_{\mathrm{sym}}/3$, as
expected.

Similarly, injecting Eq. (9) into (LABEL:eq:kumx), one obtains the following
expression for the incompressibility modulus

$\displaystyle K_{\mathrm{UM}}(x,\delta)$ $\displaystyle=$ $\displaystyle
K^{k}(\delta)+Q^{k}(\delta)x+\frac{1}{2}Z^{k}(\delta)x^{2}+o(x^{3})\,,$ (26)

with $K^{k}=K^{p}=K^{\epsilon}$ and $Q^{k}=Q^{p}$, and where the additional
coefficient in asymmetric matter reads

$\displaystyle Z^{k}(\delta)$ $\displaystyle\equiv$ $\displaystyle
Z^{k}_{\mathrm{sat}}+Z^{k}_{\mathrm{sym}}\delta^{2}\,,$ (27)

with

$\displaystyle Z^{k}_{\mathrm{sat}}$ $\displaystyle\equiv$ $\displaystyle
Z_{\mathrm{sat}}+54K_{\mathrm{sat}}+18Q_{\mathrm{sat}}\,$ (28) $\displaystyle
Z^{k}_{\mathrm{sym}}$ $\displaystyle\equiv$ $\displaystyle
Z_{\mathrm{sym}}+54K_{\mathrm{sym}}+18Q_{\mathrm{sym}}\,.$ (29)

Remark that while $K(\delta)$ controls the isoscalar and isovector dependence
of the curvature of the energy per particle in uniform matter (9), the
incompressibility (LABEL:eq:kumx) itself is driven by the parameter
$K^{k}(\delta)$ for $x=0$. The difference between $K(\delta)$ and
$K^{k}(\delta)$ reflects the contribution of the pressure, which is non-zero
as soon as the density departs from the equilibrium density
$n_{\mathrm{eq}}^{\mathrm{UM}}$, see Eq. (4). This contribution is
unavoidable, making $K^{k}(\delta)$ the true isospin dependence of the
incompressibility Piekarewicz and Centelles (2009). Fixing
$n=n_{\mathrm{sat}}$ for instance, the parameter which controls the isospin
dependence of the incompressibility is
$K^{k}_{\mathrm{sym}}=K_{\mathrm{sym}}+6L_{\mathrm{sym}}$, and not Ksym alone.
Considering $L_{\mathrm{sym}}\approx 50$ MeV and $K_{\mathrm{sym}}\approx-100$
MeV d’Etivaux  _et al._ (2019); Sagawa  _et al._ (2019), with a lower limit
provided by the unitary limit Tews  _et al._ (2017), the parameter
$K^{k}_{\mathrm{sym}}$ is even mostly controlled by $L_{\mathrm{sym}}$, and
only moderately by $K_{\mathrm{sym}}$.

In SM the incompressibility modulus can be expressed as a series expansion in
$x$ as

$\displaystyle K_{\mathrm{SM}}(x)$ $\displaystyle=$ $\displaystyle
K_{\mathrm{sat}}+(12K_{\mathrm{sat}}+Q_{\mathrm{sat}})x$ (30) $\displaystyle+$
$\displaystyle(27K_{\mathrm{sat}}+9Q_{\mathrm{sat}})x^{2}+o(x^{3})\,,$

and we introduce a new quantity,

$\displaystyle K_{\mathrm{sat},\mathrm{eq}}$ $\displaystyle\equiv$
$\displaystyle K_{\mathrm{SM}}(x=x_{\mathrm{eq}}^{\mathrm{UM}})$ (31)
$\displaystyle=$ $\displaystyle
K_{\mathrm{sat}}+(12K_{\mathrm{sat}}+Q_{\mathrm{sat}})x_{\mathrm{eq}}^{\mathrm{UM}}+o(x^{2})\,,$

which represents the incompressibility modulus of SM for the equilibrium
density $n_{\mathrm{eq}}^{\mathrm{UM}}$:

One can show that in AM the incompressibility modulus at the equilibrium
density (26) can be expressed as,

$K_{\mathrm{eq}}\equiv
K_{\mathrm{UM}}(x_{\mathrm{eq}}^{\mathrm{UM}},\delta)=K_{\mathrm{sat},\mathrm{eq}}+K_{\mathrm{sym}}^{k}\delta^{2}+o(x^{2},\delta^{4})\,.$
(32)

In Eq. (32), the isovector term
$K_{\mathrm{sym}}^{k}$=$K_{\mathrm{sym}}^{\epsilon}$ (15) depends only on
isovector empirical parameters $L_{\mathrm{sym}}$ and $K_{\mathrm{sym}}$,
while the isoscalar term only depends on isoscalar NEPs $K_{\mathrm{sat}}$ and
$Q_{\mathrm{sat}}$, provided $x_{\mathrm{eq}}^{\mathrm{UM}}$ is known
(experimentally for instance). In order to perform comparisons with
incompressibilities in nuclei $K_{A}$, it could be relevant to express the
incompressibility modulus in AM at equilibrium density
$x=x_{\mathrm{eq}}^{\mathrm{UM}}$ as,

$\displaystyle
K_{\mathrm{eq}}=K_{\mathrm{sat}}+K_{\tau}\delta^{2}+o(x^{2},\delta^{4})\,,$
(33)

where Piekarewicz and Centelles (2009)

$K_{\tau}=K_{\mathrm{sym}}-(6+Q_{\mathrm{sat}}/K_{\mathrm{sat}})L_{\mathrm{sym}}\,.$
(34)

We choose nine Skyrme models, BSK14Goriely  _et al._ (2007), BSK16Chamel  _et
al._ (2008), F0Lesinski  _et al._ (2006), LNS5Cao  _et al._ (2006), RATPRayet
_et al._ (1982), SGIIVan Giai and Sagawa (1981), SKI2Reinhard and Flocard
(1995), SKOReinhard  _et al._ (1999), SLy5Chabanat  _et al._ (1998), which NEP
are given in Table 1. For these nine interactions, while the parameter
$K^{k}_{\mathrm{sym}}$ is positive for actual values of the NEPs, the
parameter $K_{\tau}$ controlling the isovector dependence of $K_{\mathrm{eq}}$
is negative since $Q_{\mathrm{sat}}/K_{\mathrm{sat}}\approx-1.5$ from Table 1.
Note however that the value of $Q_{\mathrm{sat}}$ has never been measured and
its actual value is not necessarily in the range given in Table 1. Aside from
the finite size contribution, the ISGMR in finite nuclei is mostly correlated
with $K_{\mathrm{eq}}$, whose isospin dependence is given by $K_{\tau}$
Piekarewicz and Centelles (2009). This is the reason why the isovector
dependence of the ISGMR across isotopic chains has been correlated with the
parameter $K_{\tau}$ Li  _et al._ (2010).

It is clear from the definition of $K_{\tau}$ (34) that here also, a precise
experimental determination of $K_{\tau}$ does not necessarily lead to a better
value for the NEP $K_{\mathrm{sym}}$, since $K_{\tau}$ is mostly correlated
with $L_{\mathrm{sym}}$, which is not precisely known. In order to extract
$K_{\mathrm{sym}}$ from experimental investigations, one has to precisely know
the values of $L_{\mathrm{sym}}$ and $Q_{\mathrm{sat}}$. Note that with about
10% accuracy Khan  _et al._ (2012); Khan and Margueron (2013); Khan  _et al._
(2013), the NEP $K_{\mathrm{sat}}$ is sufficiently well known in the present
case.

An illustration of the different points where the incompressibility has been
introduced is shown in Fig. 1. It displays the behavior of the equilibrium
density as a function of $\delta$, on the example of the BSk12 functional. The
role of the incompressibility at various densities and isospin is also
displayed on the figure. It shows that several incompressibilities at various
densities are probed when the GMR is measured in a given nuclei. For instance,
in addition to the saturation density, their typical mean density is around
0.11 fm-3 Khan  _et al._ (2012); Khan and Margueron (2013). It should be noted
that $K^{k}_{\mathrm{sym}}$ drives the isospin dependence of the
incompressibility, independently of the considered density, from $x=0$
(saturation point) to $x=x^{\mathrm{UM}}_{\mathrm{eq}}$ (equilibrium point).

Figure 1: Equilibrium points for the BSk12 functional (using its Lsym and Ksat
values in Eq. (22)) in the (density, isospin asymmetry) map. The corresponding
incompressibilities are schematically indicated. The upper (lower) box are
drawn around saturation (mean) densities of experimentally accessible nuclei
($\delta<$0.2).

###  II.2 An alternative representation of the nuclear matter energy and
incompressibility modulus

In this section, we explore an alternative representation of the uniform
matter properties, where the equilibrium density $n_{\mathrm{eq}}$ is taken,
in place of the saturation density $n_{\mathrm{sat}}$. This alternative
representation is equivalent to the existing one up to $\delta^{2}$, but
generates non-quadratic terms. In the view of constraining uniform matter
parameters from measurements of incompressibilities in nuclei, it may be more
relevant to consider such a representation: the equilibrium density in uniform
matter shall be closer – than the saturation density – to the average one of
the nucleus Khan  _et al._ (2012); Khan and Margueron (2013).

In this alternative approach, the associated density parameter is set to be
$\tilde{x}=(n-n_{\mathrm{eq}}^{\mathrm{UM}})/(3n_{\mathrm{eq}}^{\mathrm{UM}})$,
from which the density $n$ is obtained as
$n/n_{\mathrm{eq}}^{\mathrm{UM}}=1+3\tilde{x}$.

The alternative density parameter $\tilde{x}$ can be expressed in term of $x$
as,

$\displaystyle\tilde{x}=\frac{x+(L_{\mathrm{sym}}/K_{\mathrm{sat}})\delta^{2}}{1-3(L_{\mathrm{sym}}/K_{\mathrm{sat}})\delta^{2}}\,.$
(35)

Similarly to Eq. (9) one can expand the energy per particle in term of
$\tilde{x}$ as

$e_{\mathrm{UM}}(\tilde{x},\delta)=\tilde{E}(\delta)+\frac{1}{2}\tilde{K}(\delta)\tilde{x}^{2}+\frac{1}{6}\tilde{Q}(\delta)\tilde{x}^{3}+\frac{1}{24}\tilde{Z}(\delta)\tilde{x}^{4}+\dots\,,$
(36)

with

$\displaystyle\tilde{E}(\delta)$ $\displaystyle=$ $\displaystyle
e(\tilde{x}=0,\delta)=E_{\mathrm{sat}}+E_{\tau}\delta^{2}\,,$ (37)
$\displaystyle\tilde{K}(\delta)$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}e(\tilde{x},\delta)}{\partial\tilde{x}^{2}}\Bigr{\rvert}_{\tilde{x}=0}=K_{\mathrm{sat}}+K_{\tau}\delta^{2}=K_{\mathrm{eq}}\,,$
(38) $\displaystyle\tilde{Q}(\delta)$ $\displaystyle=$
$\displaystyle\frac{\partial^{3}e(\tilde{x},\delta)}{\partial\tilde{x}^{3}}\Bigr{\rvert}_{\tilde{x}=0}=Q_{\mathrm{sat}}+Q_{\tau}\delta^{2}\,,$
(39) $\displaystyle\tilde{Z}(\delta)$ $\displaystyle=$
$\displaystyle\frac{\partial^{4}e(\tilde{x},\delta)}{\partial\tilde{x}^{4}}\Bigr{\rvert}_{\tilde{x}=0}=Z_{\mathrm{sat}}+Z_{\tau}\delta^{2}\,.$
(40)

It should be noted that Kτ in Eq. (38) corresponds to (34), because it is the
incompressibility at the equilibrium density, namely $\tilde{K}(\delta)$=Keq.

Imposing the equality between the $\delta^{2}$ terms in the series expansions
(9) and (36) orders by orders in $x$, one obtains the following relations,

$\displaystyle E_{\tau}$ $\displaystyle=$ $\displaystyle E_{\mathrm{sym}},$
(41) $\displaystyle K_{\tau}$ $\displaystyle=$ $\displaystyle
K_{\mathrm{sym}}-L_{\mathrm{sym}}(6K_{\mathrm{sat}}+Q_{\mathrm{sat}})/K_{\mathrm{sat}},$
(42) $\displaystyle Q_{\tau}$ $\displaystyle=$ $\displaystyle
Q_{\mathrm{sym}}-L_{\mathrm{sym}}(9Q_{\mathrm{sat}}+Z_{\mathrm{sat}})/K_{\mathrm{sat}},$
(43) $\displaystyle Z_{\tau}$ $\displaystyle=$ $\displaystyle
Z_{\mathrm{sym}}-L_{\mathrm{sym}}(12Z_{\mathrm{sat}}+Y_{\mathrm{sat}})/K_{\mathrm{sat}},$
(44)

where $Y_{\mathrm{sat}}$ is the fifth order NEP. Note that these equations are
a generalization of Eq. (16) of Ref. Piekarewicz and Centelles (2009) up to
the 4th order, and hence, the above equation for Kτ is the same than the one
of the previous subsection.

Eqs. (9) and (36) are identical up to terms in $\delta^{2}$. In Eq. (36) there
are however non-quadratic terms, which are small even when $\delta\sim 1$. The
contribution of these non-quadratic terms (because of the denominator in Eq.
(35)) is even more suppressed by the fact that finite nuclei do not explore
large values for $\delta$, since $|\delta_{A}|<0.3$, as previously discussed.
So it is possible to use both Eq. (9) or Eq. (36) to describe the energy in
finite nuclei.

Expressing the incompressibility modulus in asymmetric matter (LABEL:eq:kumx)
as a function of the density parameter $\tilde{x}$:

$K_{UM}(\tilde{x},\delta)=6(1+3\tilde{x})\frac{\partial
e(\tilde{x},\delta)}{\partial\tilde{x}}+(1+3\tilde{x})^{2}\frac{\partial^{2}e(\tilde{x},\delta)}{\partial\tilde{x}^{2}},$
(45)

where we have used $(1+3x)\partial/\partial
x=(1+3\tilde{x})\partial/\partial\tilde{x}$, allows to derive the following
expression for the incompressibility:

$\displaystyle K_{UM}(\tilde{x},\delta)$ $\displaystyle=$
$\displaystyle\tilde{K}(\delta)+[12\tilde{K}(\delta)+\tilde{Q}(\delta)]\tilde{x}$
(46)
$\displaystyle+[27\tilde{K}(\delta)+9\tilde{Q}(\delta)+\frac{1}{2}\tilde{Z}(\delta)]\tilde{x}^{2}+o(\tilde{x}^{3})\,.$

Eq. (46) provides a series expansion of the incompressibility modulus in
asymmetric matter up to $\tilde{x}^{2}$ and $\delta^{2}$, which is convenient
to use when constraining the uniform matter incompressibility from
measurements in nuclei.

We will use the alternate representation developed in this section and
confront it with the standard expansion of the nuclear matter energy in Sec. V
where we present our analysis of the speed of sound.

###  II.3 Constraints on $K_{\mathrm{sym}}$

From the existence of a lower bound on the energy of NM, on the basis of
unitary-gas considerations, the following constraint on $K_{\mathrm{sym}}$ was
obtained Tews  _et al._ (2017),

$K_{\mathrm{sym}}\approx-306.0+3.41L_{\mathrm{sym}}\pm 28.3\hbox{ MeV }\,,$
(47)

when models with $K_{\mathrm{sat}}>275$ MeV are excluded. Considering
$L_{\mathrm{sym}}\approx 50$ MeV for instance, this constraint imposes
$K_{\mathrm{sym}}\geqslant\simeq-150$ MeV (see Table 1).

In Eq. (47), the coefficients of the correlation are obtained from a fit to a
given set of model realizations. In the following, we will demonstrate the
existence of a lower limit from purely theoretical considerations.

It is possible to express the equilibrium density from $K_{\mathrm{sym}}$, by
solving the mechanical stability condition $\partial
e_{\mathrm{UM}}(x,\delta)/\partial x=0$, with an expansion of the energy to
$x^{3}$ and beyond $\delta^{2}$ approximation. The physical solution of this
second order equation is

$\displaystyle
x_{\mathrm{eq},2}^{\mathrm{UM}}(\delta)=\frac{K(\delta)}{Q(\delta)}\left[-1+\sqrt{1-2\frac{L_{\mathrm{sym}}Q(\delta)}{K(\delta)^{2}}\delta^{2}}\right]\,,$
(48)

satisfying the limit $n_{\mathrm{eq},2}^{\mathrm{UM}}\rightarrow
n_{\mathrm{sat}}$ as $\delta\rightarrow 0$. Eq. (48) is well defined if
$K(\delta)^{2}\geq 2L_{\mathrm{sym}}Q(\delta)\delta^{2}$ for all values of
$\delta$ for which equilibrium density is defined, which is ranging from SM to
very asymmetric matter. There is no equilibrium density in NM, but there is
still an equilibrium very close to NM. Since Eq. (48) weakly depends on
$\delta$ for isospin asymmetries close to NM, we fix $\delta=1$ in Eq. (48)
for simplicity. We then obtain: $K_{\mathrm{sym}}\geq-
K_{\mathrm{sat}}+\sqrt{2L_{\mathrm{sym}}Q(\delta=1)}$ or
$K_{\mathrm{sym}}\leq-K_{\mathrm{sat}}-\sqrt{2L_{\mathrm{sym}}Q(\delta=1)}$.
Considering typical values for the NEPs extracted from Table 1,
$\sqrt{2L_{\mathrm{sym}}Q(\delta=1)}\sim 70-100$ MeV, so the previous
condition gives $K_{\mathrm{sym}}\gtrsim-150$ MeV or
$K_{\mathrm{sym}}\lesssim-350$ MeV. Since the second case is excluded by the
constraint on $K_{\mathrm{sym}}$ given by considerations based on the unitary-
gas Tews  _et al._ (2017), we are then left with the first condition alone:

$K_{\mathrm{sym}}\geq-
K_{\mathrm{sat}}+\sqrt{2L_{\mathrm{sym}}Q(\delta=1)}\sim-150\hbox{ MeV}\,.$
(49)

Note that using the model averaged values of $Q_{\mathrm{sat}}$ and
$Q_{\mathrm{sym}}$ Margueron  _et al._ (2018a), we have the following
condition: $Q(\delta)>0$ for all $\delta$, constraining $Q_{\mathrm{sat}}$ and
$Q_{\mathrm{sym}}$, as $Q_{\mathrm{sym}}>-Q_{\mathrm{sat}}$. However, this
relation is not always satisfied, as shown in table 1.

Other estimates of $K_{\mathrm{sym}}$ from neutron stars observations have
been suggested: from X-ray thermal emission on seven LMXBs, it was found
$K_{\mathrm{sym}}=-85^{+82}_{-70}$ MeV d’Etivaux  _et al._ (2019);from the
analysis of GW170817 it was determined that $-259<K_{\mathrm{sym}}<32$ MeV
Carson  _et al._ (2019).

Using recent FRDM mass modelMöller  _et al._ (2012) and the neutron skin of
48Ca extracted from ($p,p^{\prime}$) experiments and fixing the nuclear
incompressibiliity $K_{\mathrm{sat}}=225\pm 20$ MeV from up-to-date
experimental data of ISGMR of 208Pb, it was found that
$K_{\mathrm{sym}}=-120\pm 40$ MeV Sagawa  _et al._ (2019). The constraint
obtained from FRDM mass model leads to fix $E_{\mathrm{sym}}=32.3\pm 0.5$ MeV
and $L_{\mathrm{sym}}=53.5\pm 15$ MeV. The neutron skin experiment gives
$L_{\mathrm{sym}}=42\pm 15$ MeV.

A compilation of 16 results from independent analyses of neutron star
observational data since GW170817 lead to the following expectation
$K_{\mathrm{sym}}\approx-107\pm 88$ MeV Li  _et al._ (2021). All these data
tend to point towards negative values of $K_{\mathrm{sym}}$, with a centroid
located around $-100$ MeV. The uncertainty is difficult to estimate, but a
conservative value may be around $-100$ MeV. Note that these results are
compatible with the constraint (49) that we derived.

##  III The compressible liquid drop model with a density dependence of the
surface tension

The CLDM has been originally developed on top of the liquid drop model, where
the bulk term is a constant Myers and Swiatecki (1969). In the CLDM Myers and
Swiatecki (1969); Weiss and Cameron (1969), the bulk term is density-dependent
and the density is fixed variationally by the mechanical stability condition.
In the present approach, we suggest an extension of the CLDM by introducing a
density-dependent surface term. We show that the present eCLDM could describe
accurately both the energy of finite nuclei in their ground state as well as
the ISGMR energy.

###  III.1 Density dependent surface tension

The novelty of the present work is the introduction of a density dependent
surface tension, which is expressed as,

$\displaystyle\sigma_{\mathrm{surf}}(n_{\mathrm{cl}},I_{\mathrm{cl}})=\sigma_{\mathrm{surf}}(I_{\mathrm{cl}})\left[1+a_{\mathrm{surf}}f(A_{\mathrm{cl}})x_{\mathrm{cl}}^{2}\right]\,,$
(50)

where $x_{\mathrm{cl}}=(n_{\mathrm{cl}}-n_{\mathrm{sat}})/3n_{\mathrm{sat}}$
and the parameter $a_{\mathrm{surf}}$ controls the density dependence of the
surface energy. In practice, it encodes the deviation from $n_{\mathrm{sat}}$.
It is then larger for nuclei for which $n_{\mathrm{cl}}$ is farther from
$n_{\mathrm{sat}}$, i.e. for light and intermediate mass nuclei as well as for
exotic nuclei. In Sec. A.3 we suggest a way to estimate $a_{\mathrm{surf}}$ by
using a single microscopic calculation of $K_{A}$ in 100Sn.

In Eq. (50), the function $f(A_{\mathrm{cl}})$ is defined as,

$f(A_{\mathrm{cl}})=\frac{1}{1+\exp[-(A_{\mathrm{cl}}-A_{0})/A_{w}]},$ (51)

where $A_{\mathrm{cl}}$ is the mass number of the considered nucleus. This
function has been introduced to suppress the density dependence of the surface
tension in light nuclei, where it appears to be unrealistically too large.
From a qualitative study, we suggest the following values for the parameters
of the function $f$: $A_{0}=70$ and $A_{w}=10$.

Figure 2: Surface tension for 120Sn. Continuous gray line show the result
without density dependence on $\sigma_{\mathrm{surf}}$. Dashed line shows the
result with density-dependent surface tension given by Eq. (50).

Figure 2 shows a comparison of the surface tension $\sigma_{\mathrm{surf}}$
from the CLDM and eCLDM as function of the cluster density $n_{\mathrm{cl}}$.
The figure shows a bell shape for the eCLDM due to its quadratic dependence on
$x_{\mathrm{cl}}$, in contrast with the horizontal line of the CLDM which does
not depend on the density. The eCLDM simulates a decrease of the surface
tension by about 30% from $n_{\mathrm{cl}}\approx 0.16$ fm-3 down to
$n_{\mathrm{cl}}\approx 0.10$ fm-3. So for typical values of the cluster
density, eCLDM reduces the surface tension to a large amount.

Note that we have also investigated other functionals of the density. For
instance, we have studied a correction term similar to Eq. (50), replacing
$x_{\mathrm{cl}}$ by $n_{\mathrm{cl}}$ and fixing $a_{\mathrm{surf}}$ as
described in Sec. A.3. We found that this correction changes the pressure to a
large amount, shifting up the cluster density $n_{\mathrm{cl}}$ to unrealistic
values (above $0.25$ fm-3 in some cases).

It should be noted that we have chosen the exponent of the density-dependent
term in Eq. (50) to be two. The reason is twofold: first, it approximately
satisfies the stationarity of the surface tension w.r.t the density, see Ref.
Blaizot (1980) for more details, and second, with such a power, it directly
contributes to the incompressibility modulus in finite nuclei. Note that a
correction proportional to $x_{\mathrm{cl}}$ has been suggested in Ref. Iida
and Oyamatsu (2004), and analysed in view of its impact on the neutron skin.
However, such a term does not satisfies the requested stationary of the
surface tension and does not contribute to the incompressibility in finite
nuclei.

###  III.2 Incompressibility in finite nuclei: $K_{A}$

The incompressibility $K_{A}$ in finite nuclei is defined as,

$\displaystyle K_{A}$ $\displaystyle\equiv$ $\displaystyle
9n_{\mathrm{cl}}\frac{\partial^{2}\epsilon_{A}}{\partial
n_{\mathrm{cl}}^{2}}\Bigr{|}_{A},$ (52)

with the energy density given by $\epsilon_{A}=e_{A}n_{\mathrm{cl}}$.

According to Eq. (52) by deriving twice the energy density w.r.t the cluster
density, we obtain the incompressibility in a nucleus as,

$\displaystyle K_{A}$ $\displaystyle=$ $\displaystyle
K_{\mathrm{sat}}+K_{\tau}\delta^{2}+\mathcal{C}_{\mathrm{Coul}}\frac{3}{5}\frac{e^{2}}{r_{0}}\left(8+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)Z^{2}A^{-4/3}$
(53) $\displaystyle+$ $\displaystyle\mathcal{C}_{\mathrm{surf}}\Bigg{[}8\pi
r^{2}_{\mathrm{cl}}\sigma_{\mathrm{surf}}\left(11+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)-12\pi
n_{\mathrm{cl}}r_{\mathrm{cl}}^{2}\frac{\partial\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}}\left(10+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)+36\pi
n_{\mathrm{cl}}^{2}r_{\mathrm{cl}}^{2}\frac{\partial^{2}\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}^{2}}\Biggr{]}A^{-1/3}\,.$

where $\mathcal{C}_{\mathrm{Coul}}$ and $\mathcal{C}_{\mathrm{surf}}$ are
coefficients (close to 1) optimized in order to reproduce nuclear experimental
masses. A detailed derivation of $K_{A}$ is given in App. A. The values of the
parameters used in the present work are given in Tabel 6. We can identify in
the above expression the incompressibility modulus $K_{\mathrm{sat}}$, the
isospin term $K_{\tau}$, the Coulomb and surface terms respectively. We have
arranged this expression to be comparable with Eq. (6.3) of Blaizot Blaizot
(1980). Note that the terms in the surface contribution which are proportional
to the derivative of the surface tension w.r.t. the cluster density are absent
in usual CLDM while in the eCLDM, these terms become proportional to the
constant $a_{\mathrm{surf}}$ introduced in Eq. (50).

###  III.3 Definition of the parameter $a_{\mathrm{surf}}$ and
incompressibility predictions within the eCLDM

The new parameter $a_{\mathrm{surf}}$ controlling the density dependence of
the surface tension, is fixed to reproduce the microscopic prediction for the
incompressibility $K_{A}$ in the doubly magic $N=Z$ nucleus 100Sn. The values
$a_{\mathrm{surf}}$ and the microscopic prediction from constrained Hartree-
Fock-Bogoliubov (CHFB), $K_{A,CHFB}(^{100}\hbox{Sn})$, are shown in Table 2
for the nine Skyrme interactions. The accuracy with which the microscopic
prediction is reproduced by the eCLDM is fixed to $<1$ MeV.

|  $K_{A,CHFB}(^{100}\hbox{Sn})$ |  $a_{\mathrm{surf}}$
---|---|---
|  MeV |
BSK14Goriely  _et al._ (2007) |  153.6 |  -19.95
BSK16Chamel  _et al._ (2008) |  154.4 |  -20.00
F0Lesinski  _et al._ (2006) |  142.3 |  -19.90
LNS5Cao  _et al._ (2006) |  150.7 |  -20.95
RATPRayet  _et al._ (1982) |  147.9 |  -20.85
SGIIVan Giai and Sagawa (1981) |  133.2 |  -19.55
SKI2Reinhard and Flocard (1995) |  155.2 |  -20.00
SKOReinhard  _et al._ (1999) |  139.3 |  -19.55
SLy5Chabanat  _et al._ (1998) |  142.8 |  -20.05
Table 2: For a set of Skyrme interactions, microscopic Constrained Hartree-
Fock-Bogoliubov predictions for $K_{A}$ in 100Sn used in the calibration of
the parameter $a_{\mathrm{surf}}$.

Since the parameter $a_{\mathrm{surf}}$ is found to be very stable and close
to $\sim-20$, the fit of the eCLDM is made into two steps: First, the values
of the coefficients $\mathcal{C}_{\mathrm{surf},\mathrm{sat}}$,
$\mathcal{C}_{\mathrm{surf},\mathrm{sym}}$, and $\mathcal{C}_{\mathrm{Coul}}$
are fitted to better reproduce the experimental nuclear masses, using an
initial value $a_{\mathrm{surf}}=-20$ (see Table 6), then in a second step,
the value of $a_{\mathrm{surf}}$ is accurately fixed by fitting
$K_{A,CHFB}(^{100}\hbox{Sn})$ for each of the Skyrme model (see Table 2). For
details about the microscopic CHFB approach, we refer for instance to Ref.
Khan  _et al._ (2013).

Figure 3: Comparison of $K_{A}$ for Pb an Sn isotopic chains. Continuous
(dashed) lines show results with CLDM (eCLDM). Constrained Hartree-Fock-
Bogoliubov (CHFB) calculations are shown in dots. Red triangles shows CHF
calculations for SLy5, i.e., microscopic calculations for Sn isotopes without
pairing.

In Fig. 3, we show as function of $A$, for Sn and Pb isotopes, the comparison
of the CLDM (solid lines) and eCLDM (53) (dashed lines) predictions against
the microscopic predictions (circles) for $K_{A}$, based on the constrained
Hartree-Fock-Bogoliubov (CHFB) approach, for the set of Skyrme interaction
listed in Tab. 2 (see Refs. Khan  _et al._ (2012); Khan and Margueron (2013)
for more details on the microscopic CHFB approach). We use the microscopic
radii calculated by each interaction, to transform $E_{\mathrm{ISGMR}}$ into
$K_{A}$ using Eq. (1). For SLy5, our results are identical to the ones given
in Ref. Khan  _et al._ (2013). By comparing CHF (red triangles) and CHFB (red
circles), we see that pairing contributes to reduce the shell effects around
$A\approx 140$ in Sn and makes the isotopic evolution of $K_{A}$ smoother.

As stated by Blaizot Blaizot (1980), the CLDM predictions are largely
overestimating $K_{A}$, since the surface energy is not explicitly density-
dependent. By adding the new term (50) for the surface energy (Eq. (62)), the
eCLDM reproduces the isotopic dependence of $K_{A}$ as predicted by the
microscopic CHFB approach (note that only one nucleus (100Sn) has been used
for the calibration of $a_{\mathrm{surf}}$).

In Sn isotopes, one can note a marked step for $K_{A}$ for $A\gtrsim 132$, in
microscopic predictions, which is not present in the eCLDM prediction. The
eCLDM predicts instead a continuous decrease of $K_{A}$ over the isotopic
chain. A similar feature, while not as pronounced, is observed for Pb isotopes
for $A\gtrsim 208$. Microscopically, these steps are understood as originating
from shell effects: in 132Sn, they are 12 occupied 1h11/2 states below the
Fermi level and 10 unoccupied 1h9/2 states. For 208Pb, the Fermi level is for
$1i_{13/2}$ and the next orbital is $1i_{11/2}$. The $\Delta L=0$ isoscalar
oscillation is therefore enhanced for neutron-rich systems belonging to these
two isotopic chains. Since shell effects are not present in the eCLDM, such
steps could not be described by our macroscopic approach. The decrease of
$K_{A}$ as nuclei get more and more neutron-rich is however well reproduced by
the eCLDM approach. Such a dependence on $A$ depends on the choice of the NEP,
as illustrated in the appendix section A.5.

A | Z | $e_{A}$ | $c_{A}$ | $K_{A}$ | $\bar{K}_{\mathrm{surf}}$ | $\bar{K}_{\mathrm{Coul}}$ | $E_{\mathrm{ISGMR}}$
---|---|---|---|---|---|---|---
|  | (MeV) |  | (MeV) | (MeV) | (MeV) | (MeV)
100 | 50 | -8.11 / -8.08 [-8.25] | -1.2 / 1.5 | 143.2 / 278.3 | -276.7 / 347.7 | -5.0 / -4.9 | 19.5 / 26.6
106 | 50 | -8.40 / -8.38 [-8.43] | -1.3 / 1.5 | 142.7 / 277.3 | -289.6 / 345.2 | -5.0 / -4.9 | 19.1 / 26.0
114 | 50 | -8.55 / -8.53 [-8.52] | -1.3 / 1.5 | 139.5 / 272.2 | -305.1 / 336.5 | -5.0 / -4.9 | 18.3 / 25.1 [15.9]
120 | 50 | -8.53 / -8.52 [-8.50] | -1.4 / 1.4 | 136.1 / 266.6 | -315.0 / 327.4 | -5.0 / -4.9 | 17.6 / 24.4 [15.5]
180 | 82 | -7.72 / -7.72 [-7.73] | -1.6 / 1.5 | 130.0 / 256.8 | -366.7 / 348.2 | -4.9 / -4.9 | 14.9 / 20.8
200 | 82 | -7.87 / -7.88 [-7.88] | -1.6 / 1.4 | 128.6 / 247.9 | -367.4 / 329.7 | -4.9 / -4.9 | 14.2 / 19.6
208 | 82 | -7.83 / -7.84 [-7.87] | -1.6 / 1.4 | 127.1 / 243.1 | -367.0 / 320.5 | -4.8 / -4.8 | 13.8 / 19.1 [13.5]
Table 3: Binding energies, ratio
$c_{A}=\bar{K}_{A,\mathrm{surf}}/K_{\mathrm{sat}}$, incompressibility $K_{A}$
in finite nuclei and the contributions from the surface and Coulomb terms, and
the ISGMR energies for a set of Sn and Pb isotopes, using the meta-model
version of the SLy5 nuclear interaction. On the two sides of the bar $/$ are
compared the values obtained from the eCLDM and the CLDM approaches. The
experimental values for the binding energies and for the $E_{\rm ISGMR}$, when
available, given by the 2020 Atomic Mass Evaluation (AME) table Huang  _et
al._ (2021) and by Garg et alGarg and Coló (2018), respectively, are given
inside brackets.

As a first application of the eCLDM, we compute the binding energies, the
incompressibilities $K_{A}$ and ISGMR energies for several Sn and Pb isotopes
and for the SLy5 Skyrme interaction (see Tab. 3). All these quantities are
given for both the eCLDM and the CLDM approaches. Experimental data for the
binding energies from AME2020 table Huang  _et al._ (2021) are also given. We
have also calculated the ratio $c_{A}=K_{A,\mathrm{surf}}/K_{\mathrm{sat}}$
for the set of nuclei. Interestingly, we found $c_{A}\approx-1.3$, which is
compatible with the calculations of Ref. Patra  _et al._ (2002) deduced from a
microscopic approach. Consistently with Fig. 3, the value obtained for $K_{A}$
with the eCLDM is considerably reduced, compared to the one provided from the
CLDM, illustrating the impact of the density-dependent surface energy term.
$K_{\mathrm{surf}}$ is also given in Tab. 3. The contribution of the density-
dependent surface energy term is large: it changes the sign of the term
$K_{\mathrm{surf}}$, from positive (CLDM) to negative (eCLDM). The $A$
dependence of $K_{\mathrm{surf}}$ is also strongly modified with the density-
dependent surface energy term. The value for $K_{\mathrm{Coul}}$ is not much
impacted by the density-dependent surface energy term. In addition,
$K_{\mathrm{Coul}}$ is compatible with the value extracted from the liquid
drop expansion Sagawa  _et al._ (2007) and are rather insensitive to the
nuclear interaction. Finally, we show, in the last column, the values for
ISGMR energies. Since this values are directly impacted by $K_{A}$, see Eq.
(1), the eCLDM shows a reduction for the $E_{\rm ISGMR}$ energies. Note that
this reduction makes the eCLDM results closer to the experimental values.

Tab. 3 illustrates one of the main feature of the eCLDM approach: the present
density-dependent surface energy term has a small impact on the binding
energies, but a large contribution to the incompressibility modulus $K_{A}$ in
finite nuclei. This justifies our fitting protocol previously described. It
also shows that the low order NEP could be adjusted to the nuclear mass table
quite independently to the higher order NEP which are fitted to $K_{A}$.

##  IV Confrontation to the nuclear experimental data

In this section, we confront the eCLDM to the nuclear data. To do so, we first
list the experimental data used for the analysis. By using the Markov chain
Monte Carlo (MCMC) approach, we then vary a set of NEP all together in order
to extract the best parameters set reproducing the experimental data. A
sensitivity analysis is shown in the appendix A.5 where we illustrate, in a
complementary way, the individual influence of the NEP to the prediction of
$K_{A}$.

###  IV.1 Experimental data for $K_{A}$

| $E_{\mathrm{ISGMR}}$ | $E_{\mathrm{ISGMR}}$ | $R_{A}$ | $K_{A}$
---|---|---|---|---
| (MeV) | (MeV) | (fm) | (MeV)
| from Ref. Garg and Coló (2018) | (this work) | (SLy5) | from Eq. (1)
90Zr | 17.58${}^{+0.06}_{-0.04}$ | $17.62\pm{0.07}$ | 4.256 | $135.6\pm{1.1}$
| 17.66${}^{+0.07}_{-0.07}$ |  |  |
92Zr | 17.71${}^{+0.09}_{-0.07}$ | $17.62\pm{0.12}$ | 4.293 | $138.0\pm{1.9}$
| 17.52${}^{+0.04}_{-0.04}$ |  |  |
94Zr | 15.75${}^{+0.27}_{-0.15}$ | $15.80\pm{0.21}$ | 4.330 | $112.9\pm{3.0}$
112Sn | 15.23${}^{+0.26}_{-0.14}$ | $15.69\pm{0.44}$ | 4.556 | $123.2\pm{6.9}$
| 16.10${}^{+0.10}_{-0.10}$ |  |  |
114Sn | 15.90${}^{+0.10}_{-0.10}$ | $15.90\pm{0.10}$ | $4.585$ | $128.2\pm{1.6}$
116Sn | 15.70${}^{+0.10}_{-0.10}$ | $15.70\pm{0.10}$ | 4.614 | $126.5\pm{1.6}$
118Sn | 15.60${}^{+0.10}_{-0.10}$ | $15.60\pm{0.10}$ | 4.641 | $126.4\pm{1.6}$
120Sn | 15.50${}^{+0.10}_{-0.10}$ | $15.50\pm{0.10}$ | 4.667 | $126.2\pm{1.6}$
122Sn | 15.20${}^{+0.10}_{-0.10}$ | $15.20\pm{0.10}$ | 4.691 | $122.6\pm{1.6}$
124Sn | 14.33${}^{+0.17}_{-0.14}$ | $14.72\pm{0.40}$ | 4.715 | $116.2\pm{6.3}$
| 15.10${}^{+0.10}_{-0.10}$ |  |  |
132Sn† | 14.80 | 14.80 | 4.803 | 121.8
204Pb | 13.70${}^{+0.10}_{-0.10}$ | $13.70\pm{0.10}$ | 5.516 | $137.7\pm{2.0}$
206Pb | 13.60${}^{+0.10}_{-0.10}$ | $13.60\pm{0.10}$ | 5.532 | $136.5\pm{2.0}$
208Pb | 13.50${}^{+0.10}_{-0.10}$ | $13.50\pm{0.10}$ | 5.548 | $135.3\pm{2.0}$

†Fictitious data.

Table 4: Experimental data for $E_{\mathrm{ISGMR}}$ and $K_{A}$, considered in
this work.

We aim at reproducing together the values of $K_{A}$ in 90,92Zr, 112-124Sn,
204-208Pb from Ref. Garg and Coló (2018), see Tab. 4 for detailed values. We
do not consider here the experimental GMR energy measured for 94Zr and
reported in Ref. Garg and Coló (2018), since it is very different from the one
measured in 90Zr and 92Zr. It is not possible, for our modeling to reproduce
this data, as it is shown hereafter in Fig. 7. In addition, we investigate the
role of a fictitious measurement of the GMR energy in 132Sn and explore
possible consequences for the determination of NEP.

We first report, in Tab. 4, the experimental data listed in Ref. Garg and Coló
(2018). For some nuclei there are different values obtained from different
experiments, see for instance 90Zr, 92Zr, 112Sn and 124Sn (the largest
differences between different experimental measurements are for 112Sn and
124Sn). In the following, we adopt an agnostic approach w.r.t. these data and
we then equally treat the measurements. It should be noted that we have then
re-calculated averaged centroids and standard deviations for nuclei were two
experimental values are reported, generating a new distribution summing the
individual ones. We have then determined the value for $K_{A}$ using Eq. (1),
where the total radius $R_{A}$ is provided by an CHFB calculation Bennaceur
and Dobaczewski (2005) using SLy5 Chabanat  _et al._ (1998) Skyrme
interaction. The last column in Tab. 4 gives the experimental values for
$K_{A}$ which are used in the confrontation of our eCLDM to nuclear data.

###  IV.2 Best parameter set from Markov-chain Monte Carlo approach

In this subsection, we vary a set of NEP in order to determine the best
parameters reproducing the experimental data. We first present the
experimental data and then the Markov-chain Monte Carlo (MCMC) approach we
adopt.

The confrontation between the experimental data for the incompressibility
$K_{A}$ (see tab. 4), and the model predictions, is based on the loss
functions $\chi_{K_{A}}$, which is defined as

$\chi_{K_{A}}^{2}=\frac{1}{N_{K_{A}}}\sum_{i}\left(\frac{K_{A,i}^{\exp}-K_{A,i}^{\mathrm{eCLDM}}}{\delta
K_{A,i}^{\exp}}\right)^{2},$ (54)

where $i$ runs over the following isotopes: 90,92Zr, 112-124Sn, 204-208Pb. We
also explore a fictitious data for 132Sn, since and experimental value of the
GMR centroid is currently under analysis Ota  _et al._ (2017).

| Esat | nsat | Ksat | Qsat | Zsat | Esym | Lsym | Ksym | Qsym | Zsym
---|---|---|---|---|---|---|---|---|---|---
| (MeV) | (fm-3) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV)
From Ref. Margueron  _et al._ (2018a) | -15.8$\pm$0.3 | 0.155$\pm$0.005 | 230$\pm$20 | 300$\pm$400 | -500$\pm$1000 | 32$\pm$ 2 | 60$\pm$15 | -100$\pm$100 | 0$\pm$ 400 | -500$\pm$1000
dist1f and dist2f | -15.8 | 0.155 | [210,250] | [-1800,600] | -500 | 32 | [40,60] | [-300,100] | 0 | -500
dist3f | -15.8 | 0.155 | [210,250] | [-1800,600] | -500 | 32 | [80,100] | [-300,100] | 0 | -500
Table 5: Priors for the NEP from Ref. Margueron  _et al._ (2018a) (first raw),
and priors considered in the present analysis to set-up dist1f, dist2f and
dist3f. The values given in interval $[a,b]$ imply that a flat prior is
considered in the MCMC approach. Other NEP are fixed to the indicated values.
We considered also the following parameters: $M^{*}_{\mathrm{sat}}=0.7$,
$\Delta M^{*}_{\mathrm{sat}}=-0.1$, $b_{\mathrm{sat}}=6.9$ and
$b_{\mathrm{sym}}=0$.

The eCLDM is also fine-tuned to experimental nuclear masses. The associated
loss function $\chi_{E}$ is defined as

$\chi_{E}^{2}=\frac{1}{N_{E}}\sum_{i}\left(\frac{E_{i}^{\exp}-E_{i}^{\mathrm{eCLDM}}}{\delta
E_{i}^{\exp}}\right)^{2}\,,$ (55)

where $i$ runs over a subset of experimental binding energy extracted from the
2020 AME mass table Huang  _et al._ (2021). To speed-up the computing time, we
do not consider all nuclei in the mass table, as in Ref. Grams  _et al._
(2022a) for instance, but instead we confront the mass model to a subset of
it. To do so, we picked-up one out of hundred data. We have checked that this
selection does not impact our results, as discussed below.

In the following we fix the NEP $E_{\mathrm{sat}}$, $E_{\mathrm{sym}}$, and
$n_{\mathrm{sat}}$ to their empirical expectations, as reported in Tab. 5. We
vary the other NEP, $K_{\mathrm{sat}}$, $Q_{\mathrm{sat}}$, $L_{\mathrm{sym}}$
and $K_{\mathrm{sym}}$, considering flat priors inside the boundaries given in
Tab. 5 and defining the prior loss function $\chi_{prior}$. The higher order
NEP $Q_{\mathrm{sym}}$, $Z_{\mathrm{sat}}$ and $Z_{\mathrm{sym}}$ have no
impact on the present analysis. Hence, they are fixed to values determined
from analyses of model predictions, see Ref. Margueron  _et al._ (2018b).
Their value is also given in Tab. 5. Finally, the effective mass, which is
parameterized by $M^{*}_{\mathrm{sat}}$ and $\Delta M^{*}_{\mathrm{sat}}$, is
also fixed in the present study.

The total loss function is obtained as the sum of $\chi_{K_{A}}$, $\chi_{E}$,
and $\chi_{prior}$. We explore three scenarios in the present study:

  * •

dist1 & dist1f: all known experimental data are considered for $K_{A}$
(90,92Zr, 112-124Sn and 204-208Pb) and the priors are taken flat, as given in
Tab. 5.

  * •

dist2 & dist2f: same as dist1 & dist1f but considering a fictitious value for
$K_{A}$ in 132Sn, as given in in Tab. 4.

  * •

dist3 & dist3f: same as dist2 & dist2f but considering a large prior for
$L_{\mathrm{sym}}$, as given in Tab. 5.

The difference between the cases dist$i$ and dist$i$f ($i=1$, 2, 3) are that
dist$i$f includes the fine tuning of the eCLDM to the experimental nuclear
masses while dist$i$ does not. In the following results, we observe that there
are very little differences between dist$i$ and dist$i$f, since the NEPs
($E_{\mathrm{sat}}$, $E_{\mathrm{sym}}$ and $n_{\mathrm{sat}}$) which play a
major role in the determination of the nuclear masses are not varied in the
present study.

Figure 4: Marginalized distributions for the NEP parameters in the three cases
which are considered here: dist1f (considering all experimental data), dist2f
(adding a fictitious measurement for 132Sn) and dist3f (with a prior for
$L_{\mathrm{sym}}=90\pm 10$ MeV). Note that for dist1f and dist2f the prior is
$L_{\mathrm{sym}}=50\pm 10$ MeV.

The marginalized distributions for the NEP parameters are shown in Figs. 4 and
5. The corner plot representation in Fig. 4 shows the one parameter
distributions on the diagonal and the correlation between the parameters off
the diagonal, while in Fig. 5 we show a zoom of the one parameter
distributions. We compare in the distributions obtained without the fictitious
data for 132Sn (dist1f, blue) and with this fictitious data (dist2f, red). We
also show the marginalized distribution when the slope of the symmetry energy
$L_{\mathrm{sym}}$ is taken to be large and around 90 MeV (dist3f, green), as
suggested by the analysis of PREX2 experimental data Reed  _et al._ (2021).
The Gaussian distributions in dashed lines represent the expected
distributions for these parameters from Ref. Margueron  _et al._ (2018b). They
are also given in the first raw of Tab. 5.

Let us first remark that the value for the parameter $Q_{\mathrm{sat}}$ is
very different from the expected values given in Tab. 5. The distribution for
$Q_{\mathrm{sat}}$ is very similar for the three cases: it is peaked at around
$\approx-950$ MeV with an uncertainty of about $150$-$200$ MeV. The value
extracted from an analysis of models predictions, since there are no direct
extraction from experimental data of this parameter, is expected to be quite
different: $\approx 300\pm 400$ MeV Ref. Margueron  _et al._ (2018b). These
values are extracted from an analysis over existing non-relativistic and
relativistic phenomenological approaches. However, it was already noticed that
the value of this parameter changes a lot from a type of nuclear interaction
to another: about -350 MeV in average for Skyrme models, around 0 for
relativistic mean-field (RMF) models and around 390 MeV for relativistic
Hartree-Fock (RHF) ones. There is therefore a large model dependence of
$Q_{\mathrm{sat}}$, which may be related to its correlation with
$K_{\mathrm{sat}}$ as suggested in Ref. Khan and Margueron (2013). The value
preferred by the GMR data points toward a region which is orthogonal to any
value of existing models. We can then deduce that in order to reproduce
correctly several isotopic chains from Zr to Pb, including Sn isotopes, the
required value for $Q_{\mathrm{sat}}$ is quite different from the typical
values given in phenomenological approaches. So the possible origin of the
difficulties faced by the usual phenomenological models in reproducing both
the Sn and Pb isotopes could take its origin in the values of the NEP
$Q_{\mathrm{sat}}$ in these models. To reproduce better Sn and Pb isotopes,
more flexibility shall be given to these models, in particular the breaking of
the correlation between $K_{\mathrm{sat}}$ and $Q_{\mathrm{sat}}$. For Skyrme
models, this could come with an additional density-dependent term, or the
’$t_{3}$’ kind, as suggested in Ref. Lesinski  _et al._ (2006).

Figure 5: One parameter marginalized distributions for the NEP parameters
$K_{\mathrm{sat}}$, $Q_{\mathrm{sat}}$, $L_{\mathrm{sym}}$ and
$K_{\mathrm{sym}}$. The distributions dist$i$f are shown in solid lines with
same colors as in Fig. 4. They are compared to the distributions dist$i$
(without fine tuning to experimental nuclear masses) in thin dotted lines. The
differences between dist$i$f and dist$i$ are small.

The second remark is about the role of a fictitious measurement of the GMR
energy in 132Sn. For simplicity, we assumed an accurate measurement as
$E_{\mathrm{ISGMR}}(^{132}$Sn$)=14.8$ MeV, see dist2f. An uncertainty in
$E_{\mathrm{ISGMR}}(^{132}$Sn$)$ will produce a result between the one
suggested by dist1f and dist2f, except if the measurement is lower than the
value we considered. Let us simplify the discussion of this fictitious data by
not considering such a case. The role of this fictitious data for
$E_{\mathrm{ISGMR}}(^{132}$Sn$)$ can be seen from the difference between
dist1f (blue) and dist2f (red) distributions. While the isoscalar NEP are
weakly impacted, the isovector NEP $K_{\mathrm{sym}}$ is largely impacted by
the fictitious data: such a new measurement would shift the expected value for
$K_{\mathrm{sym}}$ towards large and negative values.

Note also that the value of $L_{\mathrm{sym}}$ is not constrained by the
considered experimental values: $L_{\mathrm{sym}}$ fully explores the flat
prior without specific structure and it is also not correlated to other NEP.
There are however correlations between $K_{\mathrm{sat}}$ and
$Q_{\mathrm{sat}}$, as well as between $K_{\mathrm{sat}}$ and
$K_{\mathrm{sym}}$ and $Q_{\mathrm{sat}}$ and $K_{\mathrm{sym}}$. The
distribution for $K\mathrm{sat}$ is more peaked than the empirical expectation
(with a width of $\pm 20$ MeV). One of the reason is that there is correlation
between the parameter $a_{\mathrm{surf}}$, controlling the density dependence
of the surface energy, and $K_{\mathrm{sat}}$. In the present study we have
fixed $a_{\mathrm{surf}}=-20$, as resulting from the typical value we obtained
in the previous subsection. This parameter is however not fixed by any
experimental data and including its uncertainty may contribute to widen the
$K_{\mathrm{sat}}$ distribution. Another reason comes from the better
agreement of our model with the experimental data, in comparison to other
phenomenological approaches, e.g., Skyrme or RMF Khan  _et al._ (2013). Since
in our model we can fix the value of $Q_{\mathrm{sat}}$ independently of
$K_{\mathrm{sat}}$, it results in a better description of the experimental
$K_{A}$ values and the parameters $K_{\mathrm{sat}}$ and $Q_{\mathrm{sat}}$
are better determined, see Figs. 4 and 5. In other words, the uncertainties in
$Q_{\mathrm{sat}}$ impacts the one in $K_{\mathrm{sat}}$, as suggested in Ref.
Margueron and Gulminelli (2019). Since $Q_{\mathrm{sat}}$ is better known from
the present approach, it results that $K_{\mathrm{sat}}$ is also determined
with a better accuracy.

Figure 6: Marginalized distributions for the parameter $K_{\tau}$. Same legend
as in Figs. 4 and 5. The centroids for $K_{\tau}$ are given in the figure for
the cases dist$i$f.

We represent in Fig. 6 the marginalized distribution for the parameter
$K_{\tau}$, defined from Eq. (34) for the cases dist$i$ (think dotted lines)
and dist$i$f (thick solid lines). Without the fictitious GMR energy in 132Sn
(dist1 and dist1f) the $K_{\tau}$ distribution is quite flat, while when the
132Sn fictitious data is considered, the $K_{\tau}$ distribution is better
localized. For the value we considered including an accurate experimental
data, we obtain $K_{\tau}\approx-358\pm 40$ MeV ($K_{\tau}\approx-356\pm 50$
MeV) for $L_{\mathrm{sym}}\approx 50\pm 10$ MeV ($L_{\mathrm{sym}}\approx
90\pm 10$ MeV). Here also, we note the relative independence of the $K_{\tau}$
distribution in the parameter $L_{\mathrm{sym}}$.

Our results also differ from others if we do not consider the fictitious data
in 132Sn. The value $K_{\tau}\approx-550\pm 100$ MeV was extracted from the
analysis of the Sn isotopic chain only (from 112Sn to 124Sn Li  _et al._
(2010)). Note also the value $K_{\tau}\approx-500\pm 50$ MeV extracted from
the same experimental data, using different Skyrme Hamiltonians and RMF
Lagrangians Sagawa  _et al._ (2007). If we apply our analysis to the same data
points as in Ref. Li  _et al._ (2010); Sagawa  _et al._ (2007), then we obtain
$K_{\tau}\approx-330\pm 120$ MeV and $K_{\tau}\approx-270\pm 100$ MeV if we
impose to reproduce $K_{A}$ in Pb as well. Note however that when we consider
the fictitious data in 132Sn, the value for $K_{\tau}$ become more peaked.
This illustrates the role of isotopes with large isospin asymmetry in the
determination of $K_{\tau}$. However, for these data on exotic nuclei to be
effective, they need to be as accurate as the data obtained for stable nuclei.

Figure 7: Comparison of the best parameter set against the experimental data
for Zr, Sn and Pb isotopes.

Finally, we show in Fig. 7 the comparison between the experimental values for
$K_{A}$ and the values obtained with our best parameter set for each cases
dist$i$ (thin lines) and dist$i$f (thick lines). Our eCLDM model is able to
well reproduce the experimental points in Zr, Sn and Pb isotopes with a very
good accuracy. Once again, this is possibly due to the large negative Qsat
value, which points to a hint for solving the so-called Sn softness puzzle.
Note that the experimental point in 94Zr is out of reach from our model. The
difference between $K_{A}$ in 92Zr and 94Zr is too large to be reproduced. For
this reason, we decided not to include 94Zr in our fit. We also advocate for a
new measurement in 94Zr, since the present data is surprising.

In the case dist1 and dist1f, the best parameter sets provide a consistent
description of the experimental value in Zr, Sn and Pb isotopes. Note however
that the data in 124Sn is not very constraining in our case, since the
uncertainty is large. Therefore, the evolution of $K_{A}$ over the Sn isotopic
chain is quite flat in our model. The effect of including the fictitious data
in 132Sn with small uncertainty, forces our model to decrease $K_{A}$ as
function of $A$ in Sn isotopes (see dist2 and dist2f). The description of Pb
isotopes, while still good, is slightly deteriorated. It is however restored
with dist3 and dist3f, where a larger value for $L_{\mathrm{sym}}\approx 90\pm
10$ MeV is explored. However, these results are still exploratory and no
conclusion could be given without an accurate measurement of the GMR energy in
132Sn.

##  V Sound speed in nuclei and uniform matter

Figure 8: Correlation of the sound speed in NM at $n_{\mathrm{sat}}$ with the
parameter $K_{sym}$ for the three cases: dist1f, dist2f and dist3f. The
horizontal red band (Drischler 2021) depicts the constraint on the sound speed
from $\chi$EFT Hamiltonians calculated at $n_{\mathrm{sat}}$ Drischler  _et
al._ (2020). The vertical line at $K_{\mathrm{sym}}=-150$ MeV separates the
excluded values obtained from our Eq. (49) (left side) from the authorized
ones (right side).

The sound speed is an important property in transport models Danielewicz  _et
al._ (2002). It is interesting to address the effects of the nuclear
properties, e.g., NEPs, on the sound speed in uniform matter and in finite
nuclei. Moreover, the sound speed is an important ingredient in the
calculation of the tidal deformation in binary neutron stars Hinderer (2008);
Flanagan and Hinderer (2008); Somasundaram  _et al._ (2021b); Han and Steiner
(2019). Therefore a connection between the sound speeds in finite nuclei and
infinite matter could help in constraining NS observables from nuclear
experiments.

The sound speed $c_{s}$ in a nuclear fluid is largely determined from the
nuclear incompressibility in asymmetric matter. It is defined as Blaizot
(1980)

$c_{s}^{2}=\frac{K(n,\delta)}{9h(n,\delta)}\,,$ (56)

where $h$ is the enthalpy per particle $h=mc^{2}+e+P/n$. In uniform matter, we
determine the sound speed from the following quantities $e_{\mathrm{UM}}$ (9),
$P_{\mathrm{UM}}$ (5) and $K_{\mathrm{UM}}$ (4), while in finite nuclei, we
use $e_{A}$ (57), $P_{A}$ (65) and $K_{A}$ (71). All these quantities have
been defined in previous sections and in the appendix A.

We show in Fig. 8 the correlation between the speed-of-sound in NM at
$n_{\mathrm{sat}}$ and the NEP $K_{\mathrm{sym}}$. We have used the posterior
distributions corresponding to the three cases: dist1f, dist2f and dist3f for
all the NEP (here $K_{\mathrm{sat}}$, $Q_{\mathrm{sat}}$, $K_{\mathrm{sym}}$
and $L_{\mathrm{sym}}$). The findings of Fig. 8 suggest that a tight
constraint on the value of the sound-speed in NM at around saturation density,
could turn into a constraint of the value of $K_{\mathrm{sym}}$. With the
advent of ab-initio calculations such as $\chi$EFT Drischler  _et al._ (2020)
it is possible to determine a band for the sound-speed in NM. In Fig. 8, the
red band shows chiral EFT calculations for the sound speed in NM at
$n_{\mathrm{sat}}$ obtained in Ref. Drischler  _et al._ (2020). At
$n_{\mathrm{sat}}$, the intersection of the red band ($\chi$EFT) and blue
contour (dist1f) suggests that $-200\lesssim K_{\mathrm{sym}}\lesssim 50$ MeV.
We have performed a similar analysis at $2n_{\mathrm{sat}}$ but it does not
bring any additional information on $K_{\mathrm{sym}}$.

It should be noted, from Fig. 8, that the inclusion of a fictitious data in
132Sn (see the contour Dist2f) may contribute to the reduction of the band
width for the sound speed in NM, reducing the values for $K_{\mathrm{sym}}$ to
be $K_{\mathrm{sym}}\lesssim-100$ MeV and $c_{s}^{2}\lesssim 0.055c^{2}$
(instead of $0.06c^{2}$). The case of dist3f is even more interesting: if only
large values for $L_{\mathrm{sym}}\sim 90\pm 10$ MeV compatible with PREX2 are
authorized then the overlap between dist3f and $\chi$EFT occurs in the
forbidden region for $K_{\mathrm{sym}}$. In other word, there is no overlap
between dist3f and $\chi$EFT. The sound speed in NM therefore contributes to
exclude large values for $L_{\mathrm{sym}}$, as suggested by PREX2.

Figure 9: Similar to Fig. 8 but at half saturation density and with
$\delta=0.5$

It is also relevant to explore the correlation between the sound-speed and
$K_{\mathrm{sym}}$ in cases similar to what exists in heavy ion collisions at
the Fermi energy. In Fig. 9, we fix the density to be $n_{\mathrm{sat}}/2$ and
isospin asymmetry parameter $\delta=0.5$. Interestingly we see that, in
contrast with the NM case ($\delta=1$) shown in Fig. 8, the correlation is
negative (anti-correlation). This is due to the dominant contribution of the
pressure to the enthalpy, for which $K_{\mathrm{sym}}$ contributes to a large
extent: $K_{\mathrm{sym}}$ contributes to the first power in the density
parameter $x$ to the pressure, while only to the second power to the energy
per particle. The isoscalar contribution to the pressure is small in the
vicinity of saturation density. Since the leading order impact of
$K_{\mathrm{sym}}$ is an odd power in $x$, it has an opposite correlation
below saturation density, as compared to above. As in Fig. 8, we see that the
uncertainty in $L_{\mathrm{sym}}$ plays a large role, as can be inferred by
comparing dist3f with the other cases. The uncertainty induced by
$L_{\mathrm{sym}}$ is of similar magnitude as the one originating from
$K_{\mathrm{sym}}$. We can thus conclude that tighter constraints on both
$L_{\mathrm{sym}}$ and $K_{\mathrm{sym}}$ will reduce the uncertainty in the
sound speed.

Figure 10: The speed of sound in finite nuclei is shown (blue curve) as a
function of $A$ for Sn (top) and Pb (bottom). The black lines represent
nuclear matter results as indicated in the legend.

We now come back to finite nuclei, where FS terms also play a role in sound
speed. These FS terms impact the connection between the sound-speed in finite
nuclei and the sound-speed in nuclear matter. In Fig. 10, we show the sound-
speed in finite nuclei as a function of $A$ for two isotopic chains: Sn and
Pb. For this calculation, we have used the SLy5 interaction. In both panels,
the sound speed in infinite matter at $n_{\mathrm{sat}}$ is shown as black
horizontal lines. The solid black line represents SM. The red dashed line
represents AM with $\delta=0.2$ and for $n_{\mathrm{eq}}=0.157$ fm-3. So the
effect of asymmetry itself is to slightly increase the sound speed, while
shifting down the equilibrium density from SM to AM reduces the sound speed.
The main differences between the blue curve and the straight lines, originate
from the contribution of the FS terms in finite nuclei. There is a factor
approximately 2 between uniform matter and finite nuclei. The same difference
has been observed between $K_{\mathrm{sat}}$ and $K_{A}$, see for instance the
middle panel of Fig. 11. Interestingly, we see that the deviation between the
black and the blue lines increases with $A$ due to the fact that nuclei get
more and more neutron rich, and therefore the cluster density decreases. At
much larger $A$ (above the values shown in the figures), the FS terms finally
decrease in size and at the limit $A\rightarrow\infty$ finite nuclei and
uniform matter results do get closer.

##  VI Conclusions

In this work we have explored various ways to encode the density and isospin
asymmetry dependence of the incompressibility in nuclear matter. We have
discussed the dominant contribution of $L_{\mathrm{sym}}$ in the determination
of the isospin dependence of the incompressibility modulus. A better knowledge
of the incompressibility modulus in AM requires therefore an accurate
knowledge of $L_{\mathrm{sym}}$. In finite nuclei, by introducing an extended
CLDM (eCLDM) adding a density dependence to the surface tension proportional
to $x_{\mathrm{cl}}^{2}$, where
$x_{\mathrm{cl}}=(n_{\mathrm{cl}}-n_{\mathrm{sat}})/3n_{\mathrm{sat}}$, we
were able to provide a unified macroscopic model for nuclear masses and
incompressibility modulus. We have then rederived $K_{A}$ from the eCLDM
framework along the lines originally suggested by Blaizot Blaizot (1980). In
this way, the contribution of the new density-dependent term, in the surface
tension to $K_{A}$, is explicitly shown in the equations.

We have compared the predictions of the eCLDM for the nucleus
incompressibility $K_{A}$ with microscopic calculations and experimental data.
Thanks to the flexibility of the meta-model, a sensitivity analysis on the
impact of individual nuclear empirical parameter is made. As expected, the
isoscalar channel influences the absolute values of the energies while the
isoscalar one impact the slope of the $K_{A}$ as function of the isospin
asymmetry. A full exploration in the parameter space formed by
$K_{\mathrm{sat}}$, $Q_{\mathrm{sat}}$, $L_{\mathrm{sym}}$ and
$K_{\mathrm{sym}}$ is also performed, showing that the parameter
$Q_{\mathrm{sat}}$ must be approximately $Q_{\mathrm{sat}}\approx 950\pm 200$
MeV to reconcile the experimental GMR energies measured in Zr, SN and Pb
isotopes. Since this suggested value is different from the ones of
phenomenological forces, we then suggest a possible explanation of the origin
of the difficulties these forces faces in reproducing the experimental data on
$K_{A}$, on both Sn and Pb nuclei.

In addition, we explore the impact of a fictitious accurate measurement for
the GMR energy in 132Sn. We show that with such a measurement, the value of
$K_{\mathrm{sym}}$ and $K_{\tau}$ would be much better determined than they
are with the present data.

We have also derived two new constraints on $K_{\mathrm{sym}}$:

  * •

From the equilibrium density: $K_{\mathrm{sym}}\geq-
K_{\mathrm{sat}}+\sqrt{2L_{\mathrm{sym}}Q(\delta=1)}\sim-150$ MeV and
$Q_{\mathrm{sym}}>-Q_{\mathrm{sat}}$.

  * •

From the confrontation of our prediction for the sound speed in NM with the
$\chi$EFT, we found $-200\lesssim K_{\mathrm{sym}}\lesssim 50$ MeV. This
constraint could be more accurate if a measurement of the GMR energy in 132Sn
is known.

Let us remark that the constraint on $Q_{\mathrm{sym}}$ combined with the MCMC
exploration for $Q_{\mathrm{sat}}$ leads to the following consequence:
$Q_{\mathrm{sym}}\gtrsim 950\pm 200$ MeV.

In conclusion, the present work suggests a new way to analyze the experimental
$K_{A}$ and to extract the values of the NEP $K_{\mathrm{sat}}$,
$Q_{\mathrm{sat}}$, $L_{\mathrm{sym}}$ and $K_{\mathrm{sym}}$, which are the
most influential ones. This method is comparable to the microscopic Hartree-
Fock one, except that it does not describe shell effects. These shell effects
are however reduced by the treatment of the pairing, as shown in the
microscopic Hartree-Fock Bogoliubov calculations Garg and Coló (2018). The
advantage of our method is that we use the flexible nuclear meta-model to
simulate the role of the nuclear interaction. At variance with
phenomenological forces, the nuclear meta-model is able to freely choose the
best NEP which describe the experimental data. We found that the data favors a
large and negative value for $Q_{\mathrm{sat}}$ which is not possible with
phenomenological forces. We then suggest a possible origin for the observed
limitations of these forces.

We note that our ability to extract information on $Q_{\mathrm{sat}}$ and
$K_{\mathrm{sym}}$ from finite nuclei, is based on the variation in density
and isospin asymmetry explored by the isotopes defining the loss function
$\chi_{E}$. For simplicity we have based our analysis on the results of an
eCLDM where the densities and isospin asymmetries are taken flat in finite
nuclei. This is clearly an important feature which has to be improved in the
future. One may think for instance in implementing the meta-model in a
modeling of finite nuclei with better density profiles compared to the eCLDM.
Further works in this direction are therefore envisioned.

###### Acknowledgements.

J.M. and R.S. are supported by CNRS grant PICS-08294 VIPER (Nuclear Physics
for Violent Phenomena in the Universe) and the CNRS IEA-303083 BEOS project.
All authors are grateful to the CNRS/IN2P3 NewMAC master-project, and benefit
from PHAROS COST Action MP16214. This work is supported by the LABEX Lyon
Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its
financial support within the program Investissements d’Avenir
(ANR-11-IDEX-0007) of the French government operated by the National Research
Agency (ANR). GG is supported by Fonds de la Recherche Scientifique
(F.R.S.-FNRS) and the Fonds Wetenschappelijk Onderzoek - Vlaanderen (FWO)
under the Excellence of Science (EOS) Projects nr O022818F and O000422F.

##  Appendix A Derivation of the incompressibility $K_{A}$ in finite nuclei

In this appendix we derive $K_{A}$ in finite nuclei and obtain an expression
similar to the one obtained by Blaizot Blaizot (1980), but for the eCLDM
approach and where we have introduced the NEP explicitly. We detail the
derivation step by step, starting with the definition for the nuclear binding
energy, going to the pressure in the nucleus to have in the end a clear
expression for the incompressibility. In the last section of this appendix we
take advantage of the present approach to analyse the impact of the NEP to
reproduce $K_{A}$.

###  A.1 Energy per particle: $e_{A}$

We define the binding energy $e_{A}$ for the nucleus $A$ in the CLDM as,

$e_{A}\equiv e_{A,\mathrm{UM}}+e_{A,\mathrm{FS}}$ (57)

where the uniform matter energy $e_{A,\mathrm{UM}}$ is defined from the
symmetric matter and symmetry energy terms $e_{\mathrm{SM}}$ (6) and
$e_{\mathrm{sym}}$ (7) as
$e_{A,\mathrm{UM}}(n_{\mathrm{cl}},\delta_{\mathrm{cl}})\equiv
e_{\mathrm{SM}}(n_{\mathrm{cl}})+e_{\mathrm{sym}}(n_{\mathrm{cl}})\delta_{\mathrm{cl}}^{2}$,
where $n_{\mathrm{cl}}$ and $\delta_{\mathrm{cl}}$ are the equilibrium density
and isospin asymmetry $\delta_{\mathrm{cl}}=(N-Z)/A$ of a given nucleus. The
density $n_{\mathrm{cl}}$ is obtained assuming that the nucleus is at
mechanical equilibrium, i.e. $P_{A}=0$, see discussion in the next subsection.

The finite size contribution is defined as,

$\displaystyle e_{A,\mathrm{FS}}$ $\displaystyle\equiv$ $\displaystyle
e_{A,\mathrm{surf}}+e_{A,\mathrm{Coul}}$ (58)

where we consider only the surface and Coulomb terms in the present work. The
contributions originating from higher order terms in the leptodermous
expansion are disregarded in this analysis, where we present an eCLDM with a
density-dependent surface energy. However, they shall be studied in a future
work. Note also that by considering only the FS terms as given in Eq. (58),
our equations are consistent with the seminal paper by Blaizot Blaizot (1980)
(see for example Eq. (2.17)).

Considering the direct Coulomb contribution only, as well as a uniform charge
distribution in the nucleus, the Coulomb energy reads,

$\displaystyle e_{A,\mathrm{Coul}}$ $\displaystyle=$
$\displaystyle\mathcal{C}_{\mathrm{Coul}}\frac{3}{5}\frac{Z^{2}e^{2}}{R_{A}}\frac{1}{A},$
(59)

where the nucleus radius is $R_{A}=r_{\mathrm{cl}}A^{1/3}$ and
$r_{\mathrm{cl}}=(3/4\pi n_{\mathrm{cl}})^{1/3}$. The parameter
$\mathcal{C}_{\mathrm{Coul}}$, which is fitted on experimental nuclear masses
(see Ref. Grams  _et al._ (2022a) for details on the fit procedure),
represents an effective way to incorporate the effect of exchange as well as
of the surface, on the Coulomb energy. In the present fit, the experimental
masses are corrected by the odd-even mass staggering as $\tilde{E}_{\rm
ex}^{i}=E_{\rm ex}^{i}-\Delta E_{\rm ex}^{i}$, with

$\Delta E_{\rm
ex}^{i}=\left[\Delta_{\mathrm{sat}}+\Delta_{\mathrm{sym}}\left(\frac{N_{i}-Z_{i}}{A_{i}}\right)^{2}\right]\,A_{i}^{-1/3}\,\delta(N,Z)\,.$
(60)

where $\delta(N,Z)=1$ if $N$ and $Z$ are odd, $0$ if either $N$ or $Z$ is odd,
and $-1$ if both $N$ and $Z$ are even Bohr and Mottelson (1969). The
parameters $\Delta_{\mathrm{sat}}$ and $\Delta_{\mathrm{sym}}$ are varied
together with the CLDM parameters $\mathcal{C}_{i}$ in the fit to the
experimental masses. We show the optimal CLDM parameters $\mathcal{C}_{i}$ and
the odd-even mass staggering parameters for each Skyrme model in Table 6. Note
that the values we obtain for $\Delta_{\mathrm{sat}}$ and
$\Delta_{\mathrm{sym}}$ are similar to the ones determined in Ref. Vogel  _et
al._ (1984).

Model | BSK14 | BSK16 | F0 | LNS5 | RATP | SGII | SKI2 | SKO | SLy5
---|---|---|---|---|---|---|---|---|---
$\mathcal{C}_{\mathrm{Coul}}$ | 0.93/0.94 | 0.95/0.95 | 0.93/0.94 | 0.91/0.91 | 0.95/0.95 | 0.92/0.92 | 0.93/0.93 | 0.93/0.93 | 0.94/0.94
$\mathcal{C}_{\mathrm{surf}}$ | 1.03/1.02 | 1.07/1.06 | 1.09/1.08 | 0.98/0.97 | 1.07/1.06 | 0.97/0.96 | 1.00/1.00 | 1.03/1.03 | 1.07/1.06
$\mathcal{C}_{\mathrm{surf},\mathrm{sym}}$ | 0.98/0.94 | 0.92/0.87 | 1.30/1.24 | 0.93/0.90 | 0.81/0.76 | 0.58/0.54 | 1.40/1.45 | 1.26/1.25 | 1.31/1.25
$\Delta_{\mathrm{sat}}$ (MeV) | 12.5/12.4 | 12.1/12.0 | 12.5/12.4 | 12.8/12.7 | 11.9/11.8 | 12.1/12.0 | 13.3/13.4 | 13.0/12.9 | 12.2/12.1
$\Delta_{\mathrm{sym}}$ (MeV) | -37.5/-34.6 | -22.1/-19.5 | -38.3/-34.8 | -51.9/-49.8 | -14.8/-12.4 | -24.9/-21.7 | -73.2/-77.1 | -58.0/-57.2 | -42.4/-38.9
$\sqrt{\chi^{2}}$ (MeV) | 3.3/3.2 | 3.3/3.1 | 3.4/3.3 | 3.6/3.4 | 3.3/3.1 | 3.5/3.3 | 3.7/3.5 | 3.6/3.3 | 3.4/3.3
Table 6: Optimized finite size parameters and loss function $\sqrt{\chi^{2}}$
with eCLDM/CLDM. For eCLDM we use $a_{\mathrm{surf}}=-20.0$.

The surface energy is given by,

$\displaystyle e_{A,\mathrm{surf}}$ $\displaystyle=$
$\displaystyle\mathcal{C}_{\mathrm{surf}}\,4\pi\sigma_{\mathrm{surf}}R^{2}_{A}\frac{1}{A}.$
(61)

$\sigma_{\mathrm{surf,sat}}$ |  $\sigma_{\mathrm{surf,sym}}$ |  $p_{\mathrm{surf}}$
---|---|---
MeV fm-2 |  MeV fm-2 |
1.1 |  2.3 |  3.0
Table 7: Standard surface parameters for the CLDM considered in this work.
Note the associated value $b_{\mathrm{surf}}=29.9$ deduced from Eq.(64).

In the CLDM approach the surface tension is usually approximated by the
following formula Lattimer and Swesty (1991),

$\sigma_{\mathrm{surf}}(I_{\mathrm{cl}})\approx\sigma_{\mathrm{surf},\mathrm{sat}}\frac{2^{p_{\mathrm{surf}}+1}+b_{\mathrm{surf}}}{Y_{p,cl}^{-p_{\mathrm{surf}}}+b_{\mathrm{surf}}+(1-Y_{p,cl})^{-p_{\mathrm{surf}}}},$
(62)

where
$Y_{p,\mathrm{cl}}=Z_{\mathrm{cl}}/A_{\mathrm{cl}}=(1-I_{\mathrm{cl}})/2$,
$I_{\mathrm{cl}}=(N_{\mathrm{cl}}-Z_{\mathrm{cl}})/A_{\mathrm{cl}}$ and
$\sigma_{\mathrm{surf},\mathrm{sat}}$ is a parameter that determines the
surface tension of symmetric nuclei. The isospin dependence is controlled by
the parameters $b_{\mathrm{surf}}$ and $p_{\mathrm{surf}}$. Fixing the
parameter $\sigma_{\mathrm{surf},\mathrm{sat}}$ to an average value, see Tab.
7, the parameters $\mathcal{C}_{\mathrm{surf}}$ and $b_{\mathrm{surf}}$ are
fitted from the nuclear chart, while the parameter $p_{\mathrm{surf}}$ is
usually fixed to a value close to $\sim 3$ Lattimer and Swesty (1991), since
it controls the isospin dependence of the surface energy for large
asymmetries, which are not reached in finite nuclei.

For small asymmetries we could expand
$\sigma_{\mathrm{surf}}(I_{\mathrm{cl}})$ as,

$\sigma_{\mathrm{surf}}(I_{\mathrm{cl}})\approx\sigma_{\mathrm{surf},\mathrm{sat}}-\sigma_{\mathrm{surf},\mathrm{sym}}I_{\mathrm{cl}}^{2}$
(63)

with

$\sigma_{\mathrm{surf},\mathrm{sym}}=\sigma_{\mathrm{surf},\mathrm{sat}}\frac{2^{p_{\mathrm{surf}}}p_{\mathrm{surf}}(p_{\mathrm{surf}}+1)}{2^{p_{\mathrm{surf}}+1}+b_{\mathrm{surf}}}\,.$
(64)

Eq. (64) relates the parameter $b_{\mathrm{surf}}$ to the surface symmetry
energy $\sigma_{\mathrm{surf},\mathrm{sym}}$, see also Ref. Grams  _et al._
(2022a) for more details. We fit the isoscalar and isovector surface
parameters from the experimental nuclear masses. The standard surface
parameters in the CLDM approach are given in Table 7. The optimized parameters
$\mathcal{C}_{\mathrm{Coul}}$, $\mathcal{C}_{\mathrm{surf}}$ and
$\mathcal{C}_{\mathrm{surf},\mathrm{sym}}$ are given in Table 6 for the
different NEP used in the present work, together with the respective
$\sqrt{\chi^{2}}$, where $\chi^{2}=\frac{1}{N}\sum_{i=1}^{N}(E_{\rm
exp}^{i}-E_{\rm A}^{i})^{2}$. $E^{i}_{\rm exp}$ are the experimental masses,
$E^{i}_{A}$ are the predictions for the CLDM/eCLDM models for given nucleus
$i$ and $N=3375$ is the number of considered nuclei from the the 2020 Atomic
Mass Evaluation (AME) Huang  _et al._ (2021).

The novelty of the present work is the introduction of a density dependent
surface tension, see Eq. (50). It should be noted that we have chosen the
exponent of the density-dependent term in Eq. (50) to be two. The reason is
twofold: first, it approximately satisfies the stationarity of the surface
tension w.r.t the density, see Ref. Blaizot (1980) for more details, and
second, with such a power, it directly contributes to the incompressibility
modulus in finite nuclei. Note that a correction proportional to
$x_{\mathrm{cl}}$ has been suggested in Ref. Iida and Oyamatsu (2004), and
analysed in view of its impact on the neutron skin. However, such a term does
not satisfies the requested stationary of the surface tension and does not
contribute to the incompressibility in finite nuclei.

###  A.2 Pressure in finite nuclei: $P_{A}$

The pressure $P_{A}$ in finite nuclei is defined as,

$P_{A}\equiv n_{\mathrm{cl}}^{2}\frac{\partial e_{A}}{\partial
n_{\mathrm{cl}}}\Bigr{\rvert}_{A},$ (65)

which can be decomposed into a bulk term, originating from uniform matter and
a finite size contribution: $P_{A}=P_{A,\mathrm{UM}}+P_{A,\mathrm{FS}}$. The
bulk term is decomposed into SM and isospin asymmetry contributions, as in Eq.
(5): $P_{A,\mathrm{UM}}=P_{A,\mathrm{SM}}+P_{A,\mathrm{sym}}\delta^{2}$,
taking $P_{A,\mathrm{SM}}=P_{\mathrm{SM}}(n_{\mathrm{cl}})$ and
$P_{A,\mathrm{sym}}=P_{\mathrm{sym}}(n_{\mathrm{cl}})$.

Note that the functions of $R_{A}$ in the binding energy also contribute to
the pressure as,

$P_{A}=-\frac{R_{A}n_{\mathrm{cl}}}{3}\frac{\partial e_{A}}{\partial
R_{A}}\Bigr{|}_{A}\,,$ (66)

since the partial derivative w.r.t. $n_{\mathrm{cl}}$ at fixed $A$ is
equivalent to a partial derivative w.r.t. $R_{A}$, with appropriate factor,
see appendix D.

The finite size pressure is
$P_{A,\mathrm{FS}}=P_{A,\mathrm{surf}}+P_{A,\mathrm{Coul}}$, with the Coulomb
pressure term derived as

$P_{A,\mathrm{Coul}}=\frac{\mathcal{C}_{\mathrm{Coul}}}{5}\frac{Z^{2}e^{2}n_{\mathrm{cl}}}{R_{A}}\frac{1}{A}.$
(67)

The surface term is decomposed into two contributions

$P_{A,\mathrm{surf}}=P_{A,\mathrm{surf}}^{\mathrm{CLDM}}+P_{A,\mathrm{surf}}^{\mathrm{DD}}$
(68)

where the first term is the usual CLDM contribution, while the second term
originates from the new density-dependent (DD) term. They are defined as,

$\displaystyle P_{A,\mathrm{surf}}^{\mathrm{CLDM}}$ $\displaystyle=$
$\displaystyle-\frac{\mathcal{C}_{\mathrm{surf}}}{3}\,8\pi\sigma_{\mathrm{surf}}R_{A}^{2}n_{\mathrm{cl}}\frac{1}{A}\,,$
(69) $\displaystyle P_{A,\mathrm{surf}}^{\mathrm{DD}}$ $\displaystyle=$
$\displaystyle\mathcal{C}_{\mathrm{surf}}\,4\pi
R_{A}^{2}n_{\mathrm{cl}}^{2}\frac{\partial\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}}\frac{1}{A}\,.$ (70)

Note that since $\partial\sigma_{\mathrm{surf}}/\partial
n_{\mathrm{cl}}\propto x_{\mathrm{cl}}\approx 0$, the contribution of the new
DD term to the pressure is small.

Numerically, the cluster density $n_{\mathrm{cl}}$ is obtained from the
mechanical stability condition $P_{A}=0$, using the Newton-Raphson algorithm
with $n_{\mathrm{eq}}^{\mathrm{UM}}$ as starting solution.

###  A.3 Incompressibility in finite nuclei: $K_{A}$

The incompressibility $K_{A}$ in finite nuclei is defined as,

$\displaystyle K_{A}$ $\displaystyle\equiv$ $\displaystyle
9n_{\mathrm{cl}}\frac{\partial^{2}\epsilon_{A}}{\partial
n_{\mathrm{cl}}^{2}}\Bigr{|}_{A},$ (71)

with the energy density given by $\epsilon_{A}=e_{A}n_{\mathrm{cl}}$.
Similarly to the energy and the pressure, the linearity of the derivative
operator allows to decompose the incompressibility $K_{A}$ in finite nuclei as
bulk and FS terms,

$\displaystyle K_{A}$ $\displaystyle=$ $\displaystyle
K_{A,\mathrm{UM}}+K_{A,\mathrm{FS}}\,,$ (72)

where $K_{A,\mathrm{UM}}=K_{A,\mathrm{SM}}+K_{A,\mathrm{sym}}\delta^{2}$,
$K_{A,\mathrm{SM}}=K_{\mathrm{SM}}(n_{\mathrm{cl}})$ and
$K_{A,\mathrm{sym}}=K_{\mathrm{sym}}(n_{\mathrm{cl}})$.

The finite size contribution to the incompressibility are given as
$K_{A,\mathrm{FS}}=K_{A,\mathrm{surf}}+K_{A,\mathrm{Coul}}$, where the Coulomb
term reads,

$\displaystyle
K_{A,\mathrm{Coul}}=\mathcal{C}_{\mathrm{Coul}}\frac{12}{5}\frac{Z^{2}e^{2}}{R_{A}}\frac{1}{A}\,,$
(73)

and the surface term is expressed as

$\displaystyle K_{A,\mathrm{surf}}$ $\displaystyle=$
$\displaystyle\mathcal{C}_{\mathrm{surf}}\Big{[}-8\pi
R_{A}^{2}\sigma_{\mathrm{surf}}+24\pi
n_{\mathrm{cl}}R_{A}^{2}\frac{\partial\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}}$ (74) $\displaystyle+36\pi
n_{\mathrm{cl}}^{2}R_{A}^{2}\frac{\partial^{2}\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}^{2}}\Big{]}\frac{1}{A}\,.$

###  A.4 Re-expression of $K_{A}$

The first time a CLDM model was used to compute the incompressibility of
nuclei goes back to the seminal work of Blaizot Blaizot (1980). In order to
compare our expression for $K_{A}$ with his work, we dedicate this section to
re-write our equations and obtain the equivalent of Eq. (6.3) of Ref. Blaizot
(1980).

In finite nuclei, the density is different from the saturation density due to
the contribution of FS and isospin asymmetry terms. If the density parameter
$x_{\mathrm{cl}}$ remains small, Blaizot suggested to express $K_{A}$ as
Blaizot (1980),

$K_{A}=K_{\mathrm{sat}}+\delta^{2}\bar{K}_{A,\mathrm{sym}}+\bar{K}_{A,\mathrm{FS}}\,.$
(75)

The new terms $\bar{K}_{A,\mathrm{sym}}$ and $\bar{K}_{A,\mathrm{FS}}$
incorporate, in addition to the contribution $K_{A,\mathrm{sym}}$ and
$K_{A,FS}$, the shift in density between $n_{\mathrm{cl}}$ and
$n_{\mathrm{sat}}$, see Appendix B, and more specifically Eq. (91). To do so,
we consider the expression for $K_{A,\mathrm{SM}}$ up to the linear order in
$x_{\mathrm{cl}}$ from Eq. (30), where the expression for $x_{\mathrm{cl}}$ in
terms of $P_{A,\mathrm{sym}}$ and $P_{A,\mathrm{FS}}$ from Eq. (91) is
injected:

$\displaystyle K_{A,\mathrm{SM}}$ $\displaystyle=$ $\displaystyle
K_{\mathrm{sat}}-\frac{3P_{A,\mathrm{sym}}}{n_{\mathrm{cl}}K_{\mathrm{sat}}}\delta^{2}\left(12K_{\mathrm{sat}}+Q_{\mathrm{sat}}\right)$
(76) $\displaystyle\hskip
28.45274pt-\frac{3P_{A,\mathrm{FS}}}{n_{\mathrm{cl}}K_{\mathrm{sat}}}\left(12K_{\mathrm{sat}}+Q_{\mathrm{sat}}\right)\,.$

Re-ordering the different terms into $K_{A}$ gives Eq. (75) where

$\displaystyle\bar{K}_{A,\mathrm{sym}}$ $\displaystyle=$ $\displaystyle
K_{A,\mathrm{sym}}-\frac{3P_{A,\mathrm{sym}}}{n_{\mathrm{cl}}K_{\mathrm{sat}}}\left(12K_{\mathrm{sat}}+Q_{\mathrm{sat}}\right)\,,$
(77) $\displaystyle\bar{K}_{A,\mathrm{FS}}$ $\displaystyle=$ $\displaystyle
K_{A,\mathrm{FS}}-\frac{3P_{A,\mathrm{FS}}}{n_{\mathrm{cl}}K_{\mathrm{sat}}}\left(12K_{\mathrm{sat}}+Q_{\mathrm{sat}}\right)\,.$
(78)

At order $o(x_{\mathrm{cl}})$, the term $\bar{K}_{A,\mathrm{sym}}$ can be
expressed as

$\bar{K}_{A,\mathrm{sym}}\approx
K_{\mathrm{sym}}-L_{\mathrm{sym}}\left(6+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)=K_{\tau}\,.$
(79)

The FS terms could be decomposed into the Coulomb and surface contributions.
The Coulomb term reads,

$\bar{K}_{A,\mathrm{Coul}}=-\frac{3\mathcal{C}_{\mathrm{Coul}}}{5}\frac{Z^{2}e^{2}}{R_{A}}\frac{1}{A}\left(8+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)=\bar{K}_{\mathrm{Coul}}Z^{2}A^{-4/3}\,,$
(80)

with

$\displaystyle\bar{K}_{\mathrm{Coul}}$ $\displaystyle=$
$\displaystyle-\frac{3\mathcal{C}_{\mathrm{Coul}}}{5}\frac{e^{2}}{r_{0}}\left(8+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)\,.$
(81)

The surface term reads,

$\displaystyle\bar{K}_{A,\mathrm{surf}}$ $\displaystyle=$
$\displaystyle\Big{(}\bar{K}_{\mathrm{surf}}^{\mathrm{CLDM}}+\bar{K}_{\mathrm{surf}}^{\mathrm{DD},\dot{\sigma}}+\bar{K}_{\mathrm{surf}}^{\mathrm{DD},\ddot{\sigma}}\Big{)}A^{-1/3}\,.$
(82)

with

$\displaystyle\bar{K}_{\mathrm{surf}}^{\mathrm{CLDM}}$ $\displaystyle=$
$\displaystyle\mathcal{C}_{\mathrm{surf}}8\pi
r^{2}_{\mathrm{cl}}\sigma_{\mathrm{surf}}\left(11+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)\,,$
(83) $\displaystyle\bar{K}_{\mathrm{surf}}^{\mathrm{DD},\dot{\sigma}}$
$\displaystyle=$ $\displaystyle-\mathcal{C}_{\mathrm{surf}}12\pi
n_{\mathrm{cl}}r_{\mathrm{cl}}^{2}\frac{\partial\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}}\left(10+\frac{Q_{\mathrm{sat}}}{K_{\mathrm{sat}}}\right)$
(84) $\displaystyle\bar{K}_{\mathrm{surf}}^{\mathrm{DD},\ddot{\sigma}}$
$\displaystyle=$ $\displaystyle\mathcal{C}_{\mathrm{surf}}36\pi
n_{\mathrm{cl}}^{2}r_{\mathrm{cl}}^{2}\frac{\partial^{2}\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}^{2}}\,.$ (85)

Note that the first derivative term,
$\bar{K}_{\mathrm{surf}}^{\mathrm{DD},\dot{\sigma}}$, is expected to be small
since $\partial\sigma_{\mathrm{surf}}/\partial n_{\mathrm{cl}}\propto
x_{\mathrm{cl}}\approx 0$. In order to compare the above finite size
contributions for Eq. (75) with Eq. (6.3) of Blaizot Blaizot (1980), we shown
in App. C how to write $\bar{K}_{A,\mathrm{FS}}$ in Blaizot notation.

###  A.5 Sensitivity analysis

Figure 11: Incompressibility for Sn and Pb isotopes. Black dots with error
bars show the results of experimental data of Garg et al. Garg and Coló
(2018). Red dashed lines (squares) shows the predictions from eCLDM (CHFB)
with SLy5 interaction. Different line colors (light grey to black) show
variation on isoscalar empirical parameters $E_{\mathrm{sat}}$ (top left),
$n_{\mathrm{sat}}$ (top right), $K_{\mathrm{sat}}$ (bottom left) and
$Q_{\mathrm{sat}}$ (bottom right).  Figure 12: Same as Fig. 11 but for
variation on the isovector empirical parameters $E_{\mathrm{sym}}$ (top),
$L_{\mathrm{sym}}$ (center) and $K_{\mathrm{sym}}$ (bottom).

We analyse the impact of both the isoscalar NEPs ($E_{\mathrm{sat}}$,
$n_{\mathrm{sat}}$, $K_{\mathrm{sat}}$ and $Q_{\mathrm{sat}}$) in Fig. 11, and
the isovector NEPs ($E_{\mathrm{sym}}$, $L_{\mathrm{sym}}$ and
$K_{\mathrm{sym}}$), in Fig. 12. The results obtained from the microscopic
CHFB calculation based on the Skyrme SLy5 Hamiltonian Chabanat  _et al._
(1998) are shown in red square for the two figures. The experimental data of
Tab. 4, are shown in black with their error-bars.

In Fig. 11 we show the impact of the isoscalar parameters $E_{\mathrm{sat}}$
and $n_{\mathrm{sat}}$ (top), $K_{\mathrm{sat}}$ and $Q_{\mathrm{sat}}$
(bottom). The effects of $E_{\mathrm{sat}}$ and $n_{\mathrm{sat}}$ are very
small and almost unnoticeable. However, the incompressibility modulus
$K_{\mathrm{sat}}$ largely impacts $K_{A}$ with a positive correlation: the
larger $K_{\mathrm{sat}}$ the larger $K_{A}$. The impact of $Q_{\mathrm{sat}}$
is also large but less linear: there is a crossing value for $A$ for which the
impact of $Q_{\mathrm{sat}}$ is negligible. On the left of this crossing $A$,
$Q_{\mathrm{sat}}$ is correlated with $K_{A}$ and on the right of it, it is
anti-correlated. The red dashed line represents the eCLDM results using SLy5
Skyrme force. The value for $K_{\mathrm{sat}}$ which predict $K_{A}$ above the
experimental data in Sn, predict $K_{A}$ below them in Pb. It is then
difficult to fix accurately $K_{\mathrm{sat}}$ to reproduce experimental data
in both Sn and Pb isotopes. The difficulty to reproduce Sn and Pb isotopes
within the same nuclear force is indeed well know in the literature Colo‘ and
Van Giai (2004); Garg  _et al._ (2007); Patel  _et al._ (2012). However, it is
possible to use the NEP $Q_{\mathrm{sat}}$ which impact $K_{A}$, in a
different way compared to $K_{\mathrm{sat}}$, as previously commented. To
reconcile eCLDM with nuclear data, a low value for $Q_{\mathrm{sat}}$ is
preferred.

We now analyse the impact of the isovector NEPs. In Fig. 12 we plot eCLDM
predictions assuming the Skyrme SLy5 Hamiltonian (red dashed lines), and then
as for the isoscalar NEPs, we vary the NEPs one after another. As expected,
these parameters do not impact $K_{A}$ in symmetric nuclei, and have an impact
which increases as the isospin asymmetries increase. The impact of
$E_{\mathrm{sym}}$ (top panel) is however invisible at the scale of the
figure, while $L_{\mathrm{sym}}$ (middle panel) and $K_{\mathrm{sym}}$ (bottom
panel) have larger impacts: $L_{\mathrm{sym}}$ is anti-correlated with
$K_{A}$, while $K_{\mathrm{sym}}$ is correlated with $K_{A}$. As the isospin
asymmetry of the isotopes increases, the impact of these isovector NEPs gets
larger and larger. The larger decrease of $K_{A}$ as function of $A$ in Sn
isotopes, is obtained for large values of $L_{\mathrm{sym}}$, but a low value
of $K_{\mathrm{sym}}$ could also simulate the same effect.

None of the variations around the SLy5 Skyrme force seems to be preferred by
the data. It is then difficult, from this sensitivity analysis, to detect
which parameter set best reproduces the experimental nuclear data: the role of
the different NEP is complex and the values which suggest a better description
of the data, seem far from the SLy5 ones. In order to search for the best
parameter set, it is then necessary to have a more global approach, where all
the NEPS could be varied together, which is what we present in Sec. IV.2.

##  Appendix B Expression for $x_{\mathrm{cl}}$ in finite nuclei

We follow the approach of Blaizot Blaizot (1980) and rewrite $x_{\mathrm{cl}}$
as follows. From the definition of the compressibility $\chi(n)$,

$\chi=\frac{1}{n}\left(\frac{dP}{dn}\right)^{-1}\,,\;\;\;\hbox{ we have
}\;\;\;\frac{dP}{dn}=\frac{1}{n\chi}.$ (86)

In $N=Z$ nuclei, $P=P_{A,SM}$, and by integrating (86) from saturation
($n_{\mathrm{sat}}$) to equilibrium ($n_{\mathrm{cl}}$),

$P_{A,SM}(n_{\mathrm{cl}})-P_{A,SM}(n_{\mathrm{sat}})=\int^{n_{\mathrm{cl}}}_{n_{\mathrm{sat}}}\frac{1}{n\chi}dn\,.$
(87)

By definition $P_{A,SM}(n_{\mathrm{sat}})=0$, and for $x_{\mathrm{cl}}$ is
close to $n_{\mathrm{sat}}$ we approximate $\chi(n)\approx\chi(n_{0})$ with
$n_{0}\in[n_{\mathrm{sat}},n_{\mathrm{cl}}]$, leading to

$P_{A,SM}(n_{\mathrm{cl}})\approx\frac{1}{\chi(n_{0})}\log\frac{n_{\mathrm{cl}}}{n_{\mathrm{sat}}}\approx\frac{1}{\chi(n_{0})}\frac{n_{\mathrm{cl}}-n_{\mathrm{sat}}}{n_{\mathrm{sat}}}\,.$
(88)

Since $K=9/(n\chi)$, we have

$\frac{1}{\chi(n_{0})}=\frac{n_{0}K_{A,SM}(n_{0})}{9}\approx\frac{n_{\mathrm{cl}}K_{\mathrm{sat}}}{9},,$
(89)

since $K_{A,\mathrm{SM}}$ is an increasing function of the density. Finally,
we obtain

$x_{\mathrm{cl}}=\frac{n_{\mathrm{cl}}-n_{\mathrm{sat}}}{3n_{\mathrm{sat}}}=\frac{3}{n_{\mathrm{cl}}K_{\mathrm{sat}}}P_{A,SM}(n_{\mathrm{cl}}).$
(90)

Eq. (90) could be interpreted as the following: there is an equivalence
between the density shift $x_{\mathrm{cl}}$ which is different from zero for
densities different from $n_{\mathrm{sat}}$, as an effect of an external
pressure $P_{A,SM}$, shifting the equilibrium density to a slightly different
one. In finite nuclei, this extra-pressure is originating from the FS and
isospin asymmetry terms, since $P_{A}(n_{\mathrm{cl}})=0$. We therefore deduce
$P_{A,SM}(n_{\mathrm{cl}})=-P_{A,\mathrm{FS}}(n_{\mathrm{cl}})-\delta^{2}P_{A,\mathrm{sym}}(n_{\mathrm{cl}})$,
and we can rewrite Eq. (90) as

$x_{\mathrm{cl}}=-\frac{3}{n_{\mathrm{cl}}K_{\mathrm{sat}}}\left(P_{A,\mathrm{FS}}(n_{\mathrm{cl}})+\delta^{2}P_{A,\mathrm{sym}}(n_{\mathrm{cl}})\right).$
(91)

##  Appendix C Contributions to the incompressibility modulus within the
Blaizot notations

In the original notations of Blaizot Blaizot (1980), the NEP where not used,
but instead the third derivative of the energy density $\epsilon$. Using the
original notations, we obtain for the Coulomb contribution,

$\displaystyle\tilde{K}_{A,\mathrm{Coul}}$ $\displaystyle=$
$\displaystyle\frac{3}{5}\frac{Z^{2}e^{2}}{AR_{A}}\left(1-\frac{27n_{\mathrm{sat}}^{2}}{K_{\mathrm{sat}}}\frac{d^{3}\epsilon}{dn^{3}}\right)\,,$
(92)

and for the surface contribution

$\displaystyle\tilde{K}_{\mathrm{surf}}^{\mathrm{CLDM}}$ $\displaystyle=$
$\displaystyle 16\pi
r^{2}_{\mathrm{cl}}\sigma_{\mathrm{surf}}\left(1+\frac{27}{2}\frac{n_{\mathrm{sat}}^{2}}{K_{\mathrm{sat}}}\frac{d^{3}\epsilon}{dn^{3}}\Bigr{|}_{n_{\mathrm{sat}}}\right)\,,$
(93) $\displaystyle\tilde{K}_{\mathrm{surf}}^{\mathrm{DD},\dot{\sigma}}$
$\displaystyle=$ $\displaystyle-12\pi
n_{\mathrm{cl}}r_{\mathrm{cl}}^{2}\frac{\partial\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}}\left(1+27\frac{n_{\mathrm{sat}}^{2}}{K_{\mathrm{sat}}}\frac{d^{3}\epsilon}{dn^{3}}\Bigr{|}_{n_{\mathrm{sat}}}\right)\,,$
(94) $\displaystyle\tilde{K}_{\mathrm{surf}}^{\mathrm{DD},\ddot{\sigma}}$
$\displaystyle=$ $\displaystyle 36\pi
n_{\mathrm{cl}}^{2}r_{\mathrm{cl}}^{2}\frac{\partial^{2}\sigma_{\mathrm{surf}}}{\partial
n_{\mathrm{cl}}^{2}}\,.$ (95)

Where the relation between the NEP and the third derivative of the energy
density can be obtained using,

$27n_{\mathrm{sat}}^{2}\frac{d^{3}\epsilon}{dn^{3}}\Bigr{|}_{n_{\mathrm{sat}}}=9K_{\mathrm{sat}}+Q_{\mathrm{sat}}.$
(96)

##  Appendix D relation between the derivatives in $n_{\mathrm{cl}}$ and the
ones in $R_{A}$ in the eCLDM

In this section, we provide the relations between derivative as function of
$n_{\mathrm{cl}}$ and as function of $R_{A}$, considering the conservation of
the mass number $A=\frac{4}{3}\pi R_{A}^{3}n_{\mathrm{cl}}$. These relations
are employed in finite nuclei, since the FS terms have an explicit dependence
on $R_{A}$ while the bulk terms depend on $n_{\mathrm{cl}}$.

We have the following relations for the first order derivatives:

$\frac{\partial}{\partial
n_{\mathrm{cl}}}\Bigr{|}_{A}=-\frac{R_{A}}{3n_{\mathrm{cl}}}\frac{\partial}{\partial
R_{A}}\Bigr{|}_{A},\;\;\;\mathrm{\&}\;\;\;\frac{\partial}{\partial
R_{A}}\Bigr{|}_{A}=-\frac{3n_{\mathrm{cl}}}{R_{A}}\frac{\partial}{\partial
n_{\mathrm{cl}}}\Bigr{|}_{A}\,,$ (97)

and for the second derivative:

$\displaystyle\frac{\partial^{2}}{\partial
R_{A}^{2}}\Bigr{|}_{A}=12\frac{n_{\mathrm{cl}}}{R_{A}^{2}}\frac{\partial}{\partial
n_{\mathrm{cl}}}\Bigr{|}_{A}+9\frac{n_{\mathrm{cl}}^{2}}{R_{A}^{2}}\frac{\partial^{2}}{\partial
n_{\mathrm{cl}}^{2}}\Bigr{|}_{A}$ (98)

## References

  * van der Woude and Harakeh (2001) A. van der Woude and M. Harakeh,  _Giant Resonances: Fundamental High-frequency Modes of Nuclear Excitation_ (Oxford University Press, 2001).
  * Rezzolla  _et al._ (2018) L. Rezzolla, P. A. M. Pizzochero, D. I. Jones, N. Rea, and I. Vidaña,  _The Physics and Astrophysics of Neutron Stars_, Vol. 457 (Springer International Publishing, 2018).
  * Margueron  _et al._ (2018a) J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025805 (2018a).
  * Blaizot (1980) J. P. Blaizot, Phys. Rep. 64, 171 (1980).
  * Garg and Coló (2018) U. Garg and G. Coló, Progress in Particle and Nuclear Physics 101, 55 (2018).
  * Huang  _et al._ (2021) W. Huang, M. Wang, F. Kondev, G. Audi, and S. Naimi, Chinese Physics C 45, 030002 (2021).
  * Li and Xie (2021) B.-A. Li and W.-J. Xie, Phys. Rev. C 104, 034610 (2021).
  * Patra  _et al._ (2002) S. K. Patra, M. Centelles, X. Viñas, and M. Del Estal, Phys. Rev. C 65, 044304 (2002).
  * Sagawa  _et al._ (2007) H. Sagawa, S. Yoshida, G. Zeng, J. Gu, and X. Zhang, Phys. Rev. C 76, 034327 (2007).
  * Avogadro and Bertulani (2013) P. Avogadro and C. A. Bertulani, Phys. Rev. C 88, 044319 (2013).
  * Goriely  _et al._ (2009) S. Goriely, S. Hilaire, M. Girod, and S. Péru, Phys. Rev. Lett. 102, 242501 (2009).
  * Lalazissis  _et al._ (2005) G. A. Lalazissis, T. Nikšić, D. Vretenar, and P. Ring, Phys. Rev. C 71, 024312 (2005).
  * Khan  _et al._ (2012) E. Khan, J. Margueron, and I. Vidaña, Phys. Rev. Lett. 109, 092501 (2012).
  * Khan and Margueron (2013) E. Khan and J. Margueron, Phys. Rev. C 88, 034319 (2013).
  * Margueron and Gulminelli (2019) J. Margueron and F. Gulminelli, Phys. Rev. C 99, 025806 (2019).
  * Myers and Swiatecki (1969) W. D. Myers and W. J. Swiatecki, Ann. Phys. 55, 395 (1969).
  * Weiss and Cameron (1969) R. A. Weiss and A. G. W. Cameron, Can. J. Phys. 47, 2171 (1969).
  * Baym  _et al._ (1971) G. Baym, H. A. Bethe, and C. Pethick, Nucl. Phys. A 175, 225 (1971).
  * Steiner (2008) A. W. Steiner, Phys. Rev. C 77, 035805 (2008).
  * Carreau  _et al._ (2019) T. Carreau, F. Gulminelli, and J. Margueron, Eur. Phys. J. A 55, 188 (2019).
  * Grams  _et al._ (2022a) G. Grams, R. Somasundaram, J. Margueron, and S. Reddy, Phys. Rev. C 105, 035806 (2022a).
  * Grams  _et al._ (2022b) G. Grams, J. Margueron, R. Somasundaram, and S. Reddy, Eur. Phys. J. A 58, 56 (2022b).
  * Somasundaram  _et al._ (2021a) R. Somasundaram, C. Drischler, I. Tews, and J. Margueron, Phys. Rev. C 103, 045803 (2021a).
  * Goriely  _et al._ (2007) S. Goriely, M. Samyn, and J. M. Pearson, Phys. Rev. C 75, 064312 (2007).
  * Chamel  _et al._ (2008) N. Chamel, S. Goriely, and J. Pearson, Nucl Phys. A 812, 72 (2008).
  * Lesinski  _et al._ (2006) T. Lesinski, K. Bennaceur, T. Duguet, and J. Meyer, Phys. Rev. C 74, 044315 (2006).
  * Cao  _et al._ (2006) L. G. Cao, U. Lombardo, C. W. Shen, and N. V. Giai, Phys. Rev. C 73, 014313 (2006).
  * Rayet  _et al._ (1982) M. Rayet, M. Arnould, F. Tondeur, and G. Paulus, Astron. Astrophys. 116, 183 (1982).
  * Van Giai and Sagawa (1981) N. Van Giai and H. Sagawa, Physics Letters B 106, 379 (1981).
  * Reinhard and Flocard (1995) P.-G. Reinhard and H. Flocard, Nuclear Physics A 584, 467 (1995).
  * Reinhard  _et al._ (1999) P.-G. Reinhard, D. J. Dean, W. Nazarewicz, J. Dobaczewski, J. A. Maruhn, and M. R. Strayer, Phys. Rev. C 60, 014316 (1999).
  * Chabanat  _et al._ (1998) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
  * Piekarewicz and Centelles (2009) J. Piekarewicz and M. Centelles, Phys. Rev. C 79, 054311 (2009).
  * Papakonstantinou  _et al._ (2013) P. Papakonstantinou, J. Margueron, F. Gulminelli, and A. R. Raduta, Phys. Rev. C 88, 045805 (2013).
  * d’Etivaux  _et al._ (2019) N. B. d’Etivaux, S. Guillot, J. Margueron, N. Webb, M. Catelan, and A. Reisenegger, The Astrophysical Journal 887, 48 (2019).
  * Sagawa  _et al._ (2019) H. Sagawa, S. Yoshida, and L.-G. Cao, AIP Conference Proceedings 2127, 020002 (2019), https://aip.scitation.org/doi/pdf/10.1063/1.5117792 .
  * Tews  _et al._ (2017) I. Tews, J. M. Lattimer, A. Ohnishi, and E. E. Kolomeitsev, The Astrophysical Journal 848, 105 (2017).
  * Li  _et al._ (2010) T. Li, U. Garg, Y. Liu, R. Marks, B. K. Nayak, P. V. Madhusudhana Rao, M. Fujiwara, H. Hashimoto, K. Nakanishi, S. Okumura, M. Yosoi, M. Ichikawa, M. Itoh, R. Matsuo, T. Terazono, M. Uchida, Y. Iwao, T. Kawabata, T. Murakami, H. Sakaguchi, S. Terashima, Y. Yasuda, J. Zenihiro, H. Akimune, K. Kawase, and M. N. Harakeh, Phys. Rev. C 81, 034309 (2010).
  * Khan  _et al._ (2013) E. Khan, N. Paar, D. Vretenar, L.-G. Cao, H. Sagawa, and G. Colò, Phys. Rev. C 87, 064311 (2013).
  * Carson  _et al._ (2019) Z. Carson, A. W. Steiner, and K. Yagi, Phys. Rev. D 99, 043010 (2019).
  * Möller  _et al._ (2012) P. Möller, W. D. Myers, H. Sagawa, and S. Yoshida, Phys. Rev. Lett. 108, 052501 (2012).
  * Li  _et al._ (2021) B.-A. Li, B.-J. Cai, W.-J. Xie, and N.-B. Zhang, Universe 7 (2021), 10.3390/universe7060182.
  * Iida and Oyamatsu (2004) K. Iida and K. Oyamatsu, Phys. Rev. C 69, 037301 (2004).
  * Bennaceur and Dobaczewski (2005) K. Bennaceur and J. Dobaczewski, Computer Physics Communications 168, 96 (2005).
  * Ota  _et al._ (2017) S. Ota _et al._ , Riken Accel. Prog. Rep. 50, 1 (2017).
  * Margueron  _et al._ (2018b) J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025805 (2018b), arXiv:1708.06894 [nucl-th] .
  * Reed  _et al._ (2021) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021).
  * Drischler  _et al._ (2020) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125, 202702 (2020), arXiv:2004.07232 [nucl-th] .
  * Danielewicz  _et al._ (2002) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002), arXiv:nucl-th/0208016 .
  * Hinderer (2008) T. Hinderer, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Flanagan and Hinderer (2008) E. E. Flanagan and T. Hinderer, Phys. Rev. D 77, 021502 (2008), arXiv:0709.1915 [astro-ph] .
  * Somasundaram  _et al._ (2021b) R. Somasundaram, I. Tews, and J. Margueron, (2021b), arXiv:2112.08157 [nucl-th] .
  * Han and Steiner (2019) S. Han and A. W. Steiner, Phys. Rev. D 99, 083014 (2019), arXiv:1810.10967 [nucl-th] .
  * Bohr and Mottelson (1969) A. Bohr and B. R. Mottelson,  _Nuclear structure vol I_ (Addison-Wesley, 1969).
  * Vogel  _et al._ (1984) P. Vogel, B. Jonson, and P. Hansen, Physics Letters B 139, 227 (1984).
  * Lattimer and Swesty (1991) J. M. Lattimer and D. Swesty, Nuclear Phys. A 535, 331 (1991).
  * Colo‘ and Van Giai (2004) G. Colo‘ and N. Van Giai, Nuclear Physics A 731, 15 (2004).
  * Garg  _et al._ (2007) U. Garg, T. Li, S. Okumura, H. Akimune, M. Fujiwara, M. Harakeh, H. Hashimoto, M. Itoh, Y. Iwao, T. Kawabata, K. Kawase, Y. Liu, R. Marks, T. Murakami, K. Nakanishi, B. Nayak, P. Madhusudhana Rao, H. Sakaguchi, Y. Terashima, M. Uchida, Y. Yasuda, M. Yosoi, and J. Zenihiro, Nuclear Physics A 788, 36 (2007).
  * Patel  _et al._ (2012) D. Patel, U. Garg, M. Fujiwara, H. Akimune, G. Berg, M. Harakeh, M. Itoh, T. Kawabata, K. Kawase, B. Nayak, T. Ohta, H. Ouchi, J. Piekarewicz, M. Uchida, H. Yoshida, and M. Yosoi, Physics Letters B 718, 447 (2012).