# Analytical mass formula and nuclear surface properties in the ETF approximation. Part II: asymmetric nuclei

F. Aymard, F. Gulminelli, J. Margueron

► To cite this version:

F. Aymard, F. Gulminelli, J. Margueron. Analytical mass formula and nuclear surface properties in the ETF approximation. Part II: asymmetric nuclei. Journal of Physics G: Nuclear and Particle Physics, 2016, 43 (4), pp.45106. 10.1088/0954-3899/43/4/045106. in2p3-01275771

HAL Id: in2p3-01275771

https://hal.in2p3.fr/in2p3-01275771

Submitted on 17 Mar 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER

# Analytical mass formula and nuclear surface properties in the ETF approximation. Part II: asymmetric nuclei

To cite this article: François Aymard et al 2016 J. Phys. G: Nucl. Part. Phys. 43 045106

View the [article online](#) for updates and enhancements.

## Related content

- Analytical mass formula and nuclear surface properties in the ETF approximation. Part I: symmetric nuclei
François Aymard, Francesca Gulminelli and Jérôme Margueron

- Properties of the nuclear medium
M Baldo and G F Burgio

- Estimating the relevance of predictions from the Skyrme-Hartree-Fock model
P-G Reinhard

## Recent citations

- Nuclear skin and the curvature of the symmetry energy
Ad. R. Raduta and F. Gulminelli

- Constraints on the nuclear equation of state from nuclear masses and radii in a Thomas-Fermi meta-modeling approach
D. Chatterjee et al

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

**IOP Astronomy** ebooks

Part of your publishing universe and your first choice for astronomy, astrophysics, solar physics and planetary science ebooks.

iopscience.org/books/aas
# Analytical mass formula and nuclear surface properties in the ETF approximation. Part II: asymmetric nuclei

François Aymard¹, Francesca Gulminelli¹ and Jérôme Margueron²

¹CNRS and ENSICAEN, UMR6534, LPC, 14050 Caen cédex, France

²Institut de Physique Nucléaire de Lyon, Université Claude Bernard Lyon 1, IN2P3-CNRS, F-69622 Villeurbanne Cedex, France

E-mail: gulminelli@lpccaen.in2p3.fr

Received 10 November 2015, revised 17 December 2015

Accepted for publication 5 January 2016

Published 17 February 2016

<figure><img src="image_4.png" /></figure>

CrossMark

## Abstract

We have recently addressed the problem of the determination of the nuclear surface energy for symmetric nuclei in the framework of the extended Thomas-Fermi (ETF) approximation using Skyrme functionals. We presently extend this formalism to the case of asymmetric nuclei and the question of the surface symmetry energy. We propose an approximate expression for the diffuseness and the surface energy. These quantities are analytically related to the parameters of the energy functional. In particular, the influence of the different equation of state parameters can be explicitly quantified. Detailed analyses of the different energy components (local/non-local, isoscalar/iso-vector, surface/curvature and higher order) are also performed. Our analytical solution of the ETF integral improves previous models and leads to a precision of better than 200 keV per nucleon in the determination of the nuclear binding energy for dripline nuclei.

Keywords: mass model, semi-classical methods, surface symmetry energy, neutron skin, extended Thomas-Fermi

(Some figures may appear in colour only in the online journal)

## 1. Introduction

The determination of a reliable analytical or quasi-analytical mass formula, covering the whole nuclear mass table, is an important issue for applications of nuclear physics and nuclear astrophysics. The best reproduction of experimentally measured masses is obtained with
phenomenological mass formulas [1–3]. However, when the predictions concern very exotic nuclei that are far from the experimental database, purely phenomenological approaches can be misleading, and extrapolations based on microscopic models are expected to be more reliable. In particular, different parametrizations of nuclear masses fitted on density functional calculations with Skyrme forces have been proposed [4–7].

In [8], analytical modelling of the neutron and proton density profiles was proposed, which reproduces with very good accuracy Hartree–Fock (HF) ground state calculations for nuclei below and even above the neutron drip-line. Knowledge of the density profiles allows accessing the whole nuclear energy functional in the framework of the extended Thomas Fermi (ETF) approach, which is based on an expansion in powers of $\hbar$ of the energy functional [6, 9–12]. Indeed, in this semi-classical theory the non-local terms in the energy density functional are entirely replaced by local gradients.

In a recent paper [13], hereafter called Part I, we used this appealing property of ETF to develop an entirely analytical mass formula for symmetric $N = Z$ nuclei, under the hypothesis that such nuclei can be described by a single density profile, and isospin breaking effects only affect the bulk incompressibility. In particular, we have established a connection between the (mass dependent) surface tension and the different parameters of the Skyrme functional, and shown that both local and local terms have to be considered, to have a predictive model of the nuclear surface.

In this paper, we present an extension of the work of Part I to the case of asymmetric nuclei, which requires the introduction of the proton and neutron density profiles as two independent degrees of freedom. In this general case, the ETF energy integral cannot be evaluated analytically. The usual approach in the literature consists in calculating the integral numerically, with density profiles which are either parametrized [6, 8, 12], or determined with a variational calculation [10, 14–16]. The limitation of such approaches is that the decomposition of the total binding energy into its different components (isoscalar, isovector, surface, curvature, etc) out of a numerical calculation is not unambiguous nor unique [17]. Moreover, a numerical calculation makes it hard to discriminate the specific influence of the different physical parameters (equation of state properties, finite range, spin orbit, etc) on quantities such as the surface symmetry energy or the neutron skin. As a consequence, correlations between observables and physical parameters require a statistical analysis based on a large set of very different models. In this way, one may hope that the obtained correlation is not spuriously induced by the specific form of the effective interaction [18]. The correlation may also depend on several physical parameters and the statistical analysis becomes quite complex [19]. Earlier approaches in the literature introduced approximations in order to keep an analytical evaluation possible [20]. These approximations however typically neglect the presence of a neutron skin, and more generally of inhomogeneities in the isospin distribution [6]. As a consequence, the results are simple and transparent, but their validity out of the stability valley should be questioned. One of the main applications of the present work concerns the production of reliable mass tables for an extensive use in astrophysical applications [21]. For this reason, we aim at expressions which stay valid approaching the driplines. In the specific application to the neutron star inner crust, even more exotic nuclei far beyond the driplines are known to be populated [22, 23]. We will not consider this situation in the present paper, because a correct treatment of nuclei beyond the dripline imposes considering the presence of both bound and unbound states which modify the density profiles and leads to the emergence of a nucleon gas. As already stated above, optimal parametrized density profiles have been proposed for this problem [8, 17, 24], but the development of systematic approximations to analytically integrate the ETF functional in the presence of a gas is a delicate issue, which will be addressed in a forthcoming paper [25].
The plan of the paper is as follows. We first demonstrate in section 2 that a large part of the isospin dependence can be accounted for if the asymmetry dependence of the saturation density is introduced in the nuclear bulk. The residual surface symmetry part is then defined in terms of the isovector density. This energy density term is not analytically integrable, meaning that approximations have to be performed. We propose in section 3 two different approximations and critically discuss their validity in comparison both to numerical integration of the ETF functional, and to complete HF calculations using the same functional (section 3.3). The first approximation, inspired by [20], consists in neglecting the neutron skin (section 3.1). Surprisingly enough, this very crude approximation leads to an overestimation of HF energies of medium-heavy and heavy nuclei of no more than 200–400 keV/nucleon even for the most extreme dripline nuclei. Such accuracy can be obtained only if both local and non-local terms in the energy functional are separately calculated, meaning that the symmetry energy does not only depend on bulk nuclear matter properties. This might be at the origin of the well-known ambiguities in the extraction of the symmetry energy from density functional calculations of finite nuclear properties [7, 26, 29, 30]. Better accuracy for neutron rich nuclei is obtained if isospin fluctuations are accounted for, and in section 3.2 the approximation is made that the surface symmetry energy density strongly peaks at the nuclear surface. Finally, the complete mass formula is calculated for different representative Skyrme functionals in section 4. The qualitative behavior of the different energy components, that is the surface, curvature and higher order terms decomposed into isovector and isoscalar parts, and local and non-local parts, is discussed. The different analytical expressions for the mass functional are explicitly demonstrated in the appendix and can be readily used with any Skyrme interaction. The paper is summarized in section 5, and conclusions are given.

## 2. Energy decomposition and density profiles

The presence of two separate good particle quantum numbers, N and Z, implies that we have to work with a two-dimensional problem. Concerning the energy functional, it is customary to split it into an isoscalar and an isovector component:

$$ \mathcal{H}[\rho, \rho_3] = \mathcal{H}^{\text{IS}}[\rho, \rho_3 = 0] + \mathcal{H}^{\text{IV}}[\rho, \rho_3] \quad (1) $$

with:

$$ \mathcal{H}^{\text{IS}}[\rho, \rho_3] = \frac{\hbar^2}{2m}\tau + C_{\text{eff}}\rho\tau + (C_0 + C_3\rho^\alpha)\rho^2 + C_{\text{fin}}(\nabla\rho)^2 + C_{\text{so}}\mathbf{J} \cdot \nabla\rho, \quad (2) $$

$$ \mathcal{H}^{\text{IV}}[\rho, \rho_3] = \mathcal{H}^{\text{IS}}[\rho, \rho_3] - \mathcal{H}^{\text{IS}}[\rho, \rho_3 = 0] + D_{\text{eff}}\rho_3\tau_3 + (D_0 + D_3\rho^\alpha)\rho_3^2 + D_{\text{fin}}(\nabla\rho_3)^2 + D_{\text{so}}\mathbf{J}_3 \cdot \nabla\rho_3, \quad (3) $$

where we have introduced the local isoscalar and isovector particle densities, kinetic densities and spin-orbit density vectors. Detailed expressions and the definition of parameters are given in Part I.

Isoscalar densities are given by the sum of the corresponding neutron and proton densities, while isovector densities (noted with the subscript ‘3’) are given by their difference. The semi-classical Wigner-Kirkwood development in $\hbar$ allows expressing all these densities in terms of the local isoscalar $\rho = \rho_n + \rho_p$ and isovector $\rho_3 = \rho_n - \rho_p$ density profiles, as well as their gradients. In equation (2), the isoscalar energy density also depends on $\rho_3$ because of the presence of the kinetic densities $\tau = \tau_n + \tau_p$, which cannot be written as a function of $\rho$ only. Therefore, to truly obtain the isoscalar part in equation (1), we have to
consider $\mathcal{H}^{\text{IS}}[\rho, \rho_3 = 0]$. The iso-vector energy density equation (3) contains therefore terms which explicitly depend on the isovector densities, but also an isovector contribution of the so-called isoscalar component $\mathcal{H}^{\text{IS}}$.

Concerning the density profiles, we choose to work with the total density $\rho(r)$ and with the proton density profile $\rho_p(r)$. Alternatively, we could as well have used $(\rho, \rho_n)$ or $(\rho_p, \rho_n)$ as independent variables, and we have checked that these different representations lead to the same level of reproduction of full HF calculations. The density profiles are parametrized as two independent Fermi functions [8]:

$$ \rho(r) = \rho_{\text{sat}} F(r); \quad \rho_p(r) = \rho_{\text{sat},p} F_p(r). \qquad (4) $$

with:

$$ F(r) = (1 + e^{(r-R)/a})^{-1}; \quad F_p(r) = (1 + e^{(r-R_p)/a_p})^{-1}. \qquad (5) $$

In equation (4), the parameters $\rho_{\text{sat}}$ and $\rho_{\text{sat},p}$ represent the saturation densities of baryons and protons of asymmetric matter [8]. These densities depend on the asymmetry $\delta = 1 - 2\rho_{\text{sat},p}/\rho_{\text{sat}}$, which represents the nucleus bulk asymmetry, according to:

$$ \rho_{\text{sat}}(\delta) = \rho_{\text{sat}}(0) \left( 1 - \frac{3L_{\text{sym}}\delta^2}{K_{\text{sat}} + K_{\text{sym}}\delta^2} \right). \qquad (6) $$

In this expression, $K_{\text{sat}} = 9\rho_{\text{sat}}^2 \partial^2(\mathcal{H}/\rho)/\partial\rho^2|_{\rho_{\text{sat}}}$ is the nuclear (symmetric) matter incompressibility, and $L_{\text{sym}} = 3\rho_{\text{sat}}\partial(\mathcal{H}_{\text{sym}}/\rho)/\partial\rho|_{\rho_{\text{sat}}}$ and $K_{\text{sym}} = 9\rho_{\text{sat}}^2 \partial^2(\mathcal{H}_{\text{sym}}/\rho)/\partial\rho^2|_{\rho_{\text{sat}}}$ are the slope and curvature of the symmetry energy at (symmetric) saturation, where we have introduced the usual definition of the symmetry energy density:

$$ \mathcal{H}_{\text{sym}} = \frac{1}{2} \rho^2 \frac{\partial^2 \mathcal{H}}{\partial \rho_3^2} \Big|_{\rho_3=0}. \qquad (7) $$

As a consequence, the radius parameters $R, R_p$ entering equation (5) also depend on the nucleus bulk asymmetry $\delta$. They are related to the particle and proton number of the nucleus by [8]:

$$ R = R_{\text{HS}} \left[ 1 - \frac{\pi^2}{3} \left( \frac{a}{R_{\text{HS}}} \right)^2 + O \left( \left( \frac{a}{R_{\text{HS}}} \right)^4 \right) \right], \qquad (8) $$

where $R_{\text{HS}} = A^{1/3}r_{\text{sat}}(\delta)$ is the equivalent homogeneous sphere radius, and $r_{\text{sat}}(\delta) = (\frac{4}{3}\pi\rho_{\text{sat}}(\delta))^{-1/3}$ is the mean radius per nucleon. A similar relation holds for the proton number, with $R \to R_p, R_{\text{HS}} \to R_{\text{HS}_p}, a \to a_p, A \to Z, r_{\text{sat}} \to r_{\text{sat},p}, \rho_{\text{sat}} \to \rho_{\text{sat},p}$. The diffusenesses $a$ and $a_p$ will be calculated in section 3 by a minimization of the surface energy. The bulk asymmetry $\delta$ differs from the global asymmetry $I = 1 - 2Z/A$ because of the competing effect of the Coulomb interaction and symmetry energy, which act in opposite directions in determining the difference between the proton and neutron radii [14, 15, 31]:

$$ \delta = \frac{I + \frac{3a_c Z^2}{8Q A^{5/3}}}{1 + \frac{9J_{\text{sym}}}{4Q A^{1/3}}}. \qquad (9) $$

In this equation, $J_{\text{sym}} = \mathcal{H}_{\text{sym}}[\rho_{\text{sat}}]/\rho_{\text{sat}}$ is the symmetry energy per nucleon at the saturation density of symmetric matter, $Q$ is the surface stiffness coefficient, and $a_c$ is the Coulomb parameter. Because of the complex interplay between Coulomb and skin effects, the bulk asymmetry $\delta$ of a globally symmetric $I = 0$ nucleus is not zero, though small for nuclei in the
<figure><img src="image_5.png" /><figcaption>Figure 1. Left side: bulk asymmetry equation (9) as a function of the global asymmetry I for nuclei within the theoretical driplines evaluated from the SLy4 [27] energy functional. The different colors correspond to different intervals in mass number: 40 ≤ A < 100 in red, 100 ≤ A < 150 in blue, 150 ≤ A < 200 in green, A ≥ 200 in gray. The function y = x is also plotted (black). Right side: β-stable nuclei (green), unstable nuclei synthesized in the laboratory [28] (red), theoretical neutron and proton driplines evaluated from the SLy4 [27] energy functional (black squares) and some iso-δ lines (blue dots) are plotted in the N, Z plane.</figcaption></figure>

nuclear chart. We have shown in [17] that accounting for the δ dependent saturation density gives a reasonably good approximation of the isospin symmetry breaking effects in I = 0 nuclei. A complete discussion on this point can be found in [26].

As a consequence, the interval of δ is slightly smaller than the interval of I over the Periodic Table. The relation between the global asymmetry and the asymmetry in the nuclear bulk is shown in the left part of figure 1. From this figure we can see that δ is a slowly increasing function of the global asymmetry I. This value increases to −0.1 < δ < 0.3 if we consider the ensemble of the heavy and medium-heavy nuclei within the driplines<fnref n="3" />. It is also observed from figure 1 that as the mass A increases, δ becomes closer to I, as expected from the analytical expression (9). The right part of figure 1 shows in the (N, Z) plane the heavy and medium-heavy measured nuclei, the theoretical neutron and proton driplines evaluated from the SLy4 energy functional (black dots), and some iso-δ lines (blue dots). We can see that the theoretical neutron dripline well matches the iso-δ line δ ≈ 0.3, which roughly corresponds to I ≈ 0.4. This means that in the following, we will be interested in approximations producing reliable formulae up to δ ≈ 0.3.

## 2.1. Bulk energy

Following the same procedure developed in Part I for symmetric nuclei, we can define the bulk energy in asymmetric matter as:

$$ 
\begin{align} 
E_b(\delta) &= \mathcal{H}_{\text{sat}}(\delta) V_{\text{HS}}(\delta) = \lambda_{\text{sat}}(\delta) A; \nonumber \\ 
\mathcal{H}_{\text{sat}}(\delta) &= \lambda_{\text{sat}}(\delta) \rho_{\text{sat}}(\delta), \tag{10}
\end{align} 
$$

where $V_{HS}(\delta) = 4/3\pi R_{HS}^3(\delta) = A/\rho_{sat}(\delta)$ is the equivalent homogeneous sphere volume and $\lambda_{sat}(\delta)$ corresponds to the chemical potential of asymmetric nuclear matter:

<footnote n="3">³ Notice that extreme values of the order $\delta \approx 0.4$ are obtained only for the lightest nuclei, not considered in this paper.</footnote>
$$ \left. \frac{\partial \mathcal{H}}{\partial \rho} \right|_{[\rho_{\text{sat}}(\delta), \rho_{\text{sat},3}(\delta)]} = \lambda_{\text{sat}}(\delta) = \left. \frac{\mathcal{H}}{\rho} \right|_{[\rho_{\text{sat}}(\delta), \rho_{\text{sat},3}(\delta)]} . \qquad (11) $$

Here, $\rho_{\text{sat},3} = \rho_{\text{sat}} - 2\rho_{\text{sat},p}$, and the total energy density $\mathcal{H}$ is given by equation (1).

## 2.2. Decomposition of the surface energy

The surface energy $E_s(\delta)$ corresponds to finite-size effects and can be decomposed into a plane surface, a curvature, and higher order terms. It is defined as the difference between the total and the bulk $E_b(\delta)$ energy,

$$ E_s(\delta) = \int d\mathbf{r} \mathcal{H}[\rho, \rho_3] - \mathcal{H}_{\text{sat}}(\delta) V_{\text{HS}}(\delta), \qquad (12) $$

where the functional $\mathcal{H}$ depends on the two densities $\rho$ and $\rho_3 = \rho - 2\rho_p$ and on the two gradients $\nabla\rho$ and $\nabla\rho_3 = \nabla\rho - 2\nabla\rho_p$. Making again the decomposition of the energy density into an isoscalar (only depending on the total density) and an isovector component (depending on $\rho$ and $\rho_3$), we get from equations (2) and (3):

$$ E_s = E_s^{\text{IS}} + E_s^{\text{IV}}, \qquad (13) $$

with

$$ E_s^{\text{IS}} = 4\pi \int_0^\infty dr \left\{ \mathcal{H}[\rho, \rho_3 = 0] - \frac{\mathcal{H}[\rho_{\text{sat}}, \rho_{\text{sat},3} = 0]}{\rho_{\text{sat}}(\delta)} \rho \right\} r^2 \qquad (14) $$

$$ E_s^{\text{IV}} = 4\pi \int_0^\infty dr \left\{ \mathcal{H}[\rho, \rho_3] - \frac{\mathcal{H}[\rho_{\text{sat}}, \rho_{\text{sat},3}]}{\rho_{\text{sat}}(\delta)} \rho \right\} r^2 - E_s^{\text{IS}}. \qquad (15) $$

It is interesting to remark that equation (12) is not the only possible definition of the surface energy in a multi-component system. Indeed in a two-component system, there are two possible definitions of the surface energy which depend on the definition of the bulk energy in the cluster [32–34]: the first one is given by equation (12) and corresponds to identifying the bulk energy of a system of $N$ neutrons and $Z$ protons to the energy of an equivalent piece of nuclear matter. The second definition $E_s \equiv E - \mu_n N - \mu_p Z + pV$ corresponds to the grand canonical thermodynamical Gibbs definition, and gives the quantity to be minimized in the variational calculation conserving proton and neutron number. Though this second definition has been often employed in the ETF literature [16, 32–34], the first one, equation (12), is the most natural definition in the present context. Indeed, using the decomposition equation (1) between isoscalar and isovector energy densities, only this definition allows recovering analytical results for the isoscalar energy, as shown in Part I. Moreover, we have shown in [17] that the best reproduction of full HF calculations is achieved considering that the bulk energy in a finite nucleus scales with the bulk asymmetry $\delta$ as in equation (12), rather than with the total asymmetry $I$, as it is implied by the Gibbs definition.

Let us first concentrate on the isoscalar surface energy. The dependence of the surface energy on the bulk asymmetry $\delta$ implies that its decomposition into an isoscalar and an isovector part is not straightforward. Indeed, although the isoscalar energy $E_s^{\text{IS}}$ does not depend on the isospin asymmetry profile $\rho_3(r)$, it does depend on the bulk isospin asymmetry $\delta$ through the isospin dependence of the saturation density $\rho_{\text{sat}}(\delta)$ appearing in the density profile $\rho$ equation (4). Moreover, in equation (14) the isoscalar bulk term which is removed depends directly on $\delta$ too, because of the equivalent volume $V_{\text{HS}} = A/\rho_{\text{sat}}(\delta)$. The quantity $E_s^{\text{IS}}$ has, therefore, an implicit dependence on isospin asymmetry $\delta$. The isoscalar surface energy
<figure><img src="image_6.png" /><figcaption>**Figure 2.** As a function of the isospin asymmetry are shown the diffuseness (left side) and isoscalar surface energy per nucleon for four isobaric chains (right side, A = 400: full lines, A = 200: dotted lines, A = 100: dashed-dotted lines, A = 50: dashed lines). Red lines: slab diffuseness $a_{\text{slab}}$ from equation (21). Blue lines: spherical diffuseness $a_{\text{sph}}$ equation (20). Green lines: quadratic diffuseness $a_{\text{HF}}$, fitted from HF density profiles in [8].</figcaption></figure>

can be calculated exactly for any nucleus of any asymmetry, with the expressions developed in Part I. In particular we can distinguish a plane surface, a curvature, and a mass independent term:

$$E_s^{\text{IS}} = E_{\text{surf}}^{\text{IS}} + E_{\text{curv}}^{\text{IS}} + E_{\text{ind}}^{\text{IS}} + O\left(\left(\frac{a(A, \delta)}{r_{\text{sat}}(\delta)}\right)^4 A^{-1/3}\right), \quad (16)$$

where the different terms are given by:

$$E_{\text{surf}}^{\text{IS}} = \left[ C_{\text{surf}}^{\text{L}} + \frac{1}{a^2(A, \delta)} C_{\text{surf}}^{\text{NL}} \right] \left( \frac{a(A, \delta)}{r_{\text{sat}}(\delta)} \right) A^{2/3} \quad (17)$$

$$E_{\text{curv}}^{\text{IS}} = \left[ C_{\text{curv}}^{\text{L}} + \frac{1}{a^2(A, \delta)} C_{\text{curv}}^{\text{NL}} \right] \left( \frac{a(A, \delta)}{r_{\text{sat}}(\delta)} \right)^2 A^{1/3} \quad (18)$$

$$E_{\text{ind}}^{\text{IS}} = \left[ C_{\text{ind}}^{\text{L}} + \frac{1}{a^2(A, \delta)} C_{\text{ind}}^{\text{NL}} \right] \left( \frac{a(A, \delta)}{r_{\text{sat}}(\delta)} \right)^3. \quad (19)$$

The local $C_i^L$ and non-local $C_i^{NL}$ functions are given in Part I.

The isoscalar diffuseness can be variationally obtained imposing $\partial E_s / \partial a = 0$, giving the following estimation $a_{\text{sph}}$ for the a parameter:

$$3C_{\text{ind}}^{\text{L}} \left(\frac{a_{\text{sph}}}{r_{\text{sat}}}\right)^4 + 2C_{\text{curv}}^{\text{L}} A^{1/3} \left(\frac{a_{\text{sph}}}{r_{\text{sat}}}\right)^3 + \left(C_{\text{surf}}^{\text{L}} A^{2/3} + \frac{1}{r_{\text{sat}}^2} C_{\text{ind}}^{\text{NL}}\right) \left(\frac{a_{\text{sph}}}{r_{\text{sat}}}\right)^2 - \frac{1}{r_{\text{sat}}^2} C_{\text{surf}}^{\text{NL}} A^{2/3} = 0. \quad (20)$$

If curvature and A-independent terms are neglected, a particularly simple solution $a = a_{\text{slab}}$ is obtained [9, 13]:
$$a_{\text{slab}} = \sqrt{\frac{c_{\text{surf}}^{\text{NL}}}{c_{\text{surf}}^{\text{L}}}} \qquad (21)$$

Though we have considered only isoscalar terms, the diffuseness $a$ does depend on the isospin asymmetry $\delta$ because of the $\delta$ dependence of the saturation density, and this is true within both approximations ($a = a_{\text{sph}}$ or $a = a_{\text{slab}}$). These results, as well as the fit from HF density profiles [8], where mass independence and quadratic behavior in $\delta$ is assumed (that is: $a = a_{\text{HF}} = C_1 + C_2\delta^2$), are shown in figure 2. Concerning the mass-dependence of equation (20) (blue lines labelled $a_{\text{sph}}$), we observe a slight spread for masses from $A = 50$ to $A = 400$, corroborating both the mass independence assumption in the HF fit $a_{\text{HF}}$ [8]: to obtain the diffuseness we can neglect the mass dependence and limit to terms $\propto A^{2/3}$ (red line, labelled $a_{\text{slab}}$). However, one can see that the dependence found from the variational equation is opposite to the one exhibited by the fit to HF results: the diffuseness decreases with $\delta$ instead of increasing. It is difficult to believe that such a huge and qualitative difference might come from the difference between ETF and HF. The discrepancy rather suggests that the variational procedure should include the isovector energy to obtain the correct behavior of the diffuseness with the isospin asymmetry. Indeed, we will see in section 3 that adding the isovector part reverses the trend. This discussion shows that, in the case of asymmetric nuclei, equation (21), which only takes into account the isoscalar terms, is not a good approximation to find the diffuseness. This is at variance with symmetric nuclei, see Part I. This statement is confirmed by the right part of the same figure, where the isoscalar surface energy per nucleon is plotted for different isobaric chains and for different prescriptions for the diffuseness. The full red and the dashed-dotted lines stand for the diffuseness $a_{\text{slab}}$ given by equation (21) and the full expression $a_{\text{sph}}$ from equation (20), respectively. There is almost no difference in the isoscalar surface energy for these two prescriptions. In addition, the observed $\delta$ dependence is extremely weak. The isoscalar surface energy evaluated with the quadratic diffuseness $a_{\text{HF}}$ [8] is represented by the dashed green line. A qualitative and quantitative difference is observed with respect to the two other curves. This indicates again that the isoscalar and isovector component of the surface energy cannot be treated separately, and the correct $\delta$ dependence of the isoscalar surface energy, as well as of the isoscalar diffuseness, requires the consideration of the total surface energy in the variational principle. It is also interesting to analyze the $\delta$ dependence of the surface symmetry energy based on the fitted quadratic diffuseness: its sign is positive, which contrasts with studies based on liquid-drop parametrizations of the nuclear mass [5, 7, 29–31, 35]. This behavior is due to our choice of definition of the surface in a two component system, as discussed at length in [17].

### 3. Approximations for the isovector energy

In this section, we focus on the residual isovector surface part $E_s^{\text{IV}}$ defined by equation (15), which cannot be written as integrals of Fermi functions as in the previous sections. Indeed, the isovector density $\rho_3$ appearing in the energy density is not a Fermi function, meaning that it cannot be analytically integrated to evaluate $E_s^{\text{IV}}$. Approximations are needed to develop an analytic expression for this part of the energy, and we will consider in the following two different approaches. At the end, we will verify the accuracy of our final formulae, in comparing the analytical expressions with HF calculations.
### 3.1. No-skin approximation

As a first approximation, we neglect all inhomogeneities in the isospin distribution in the same spirit as [20]. This simplification consists in replacing the isospin asymmetry profile $\rho_3(r)/\rho(r)$ in equation (15) by its mean value $\langle\delta\rangle$. Within this approximation, the local isovector energy only depends on the total baryonic density profile $\rho$ defined equation (4), and the non-local isovector part, involving gradients $\nabla\rho_3$, is identically zero. In other words, this approximation amounts to neglecting the non-local contribution to the isovector surface energy.

Integrating in space the equality $\rho_3(r) = \langle\delta\rangle\rho(r)$ we immediately obtain that the mean value of the isospin distribution is given by the global asymmetry of the nucleus:

$$ \langle\delta\rangle = \frac{N-Z}{A} = I. \qquad (22) $$

In particular, in this approximation, the bulk isospin asymmetry $\delta$ is equal to the global asymmetry $I$, at variance with the more elaborated relation between $\delta$ and $I$ given by equation (9). In neglecting isospin inhomogeneities, we indeed neglect both neutron skin and Coulomb effects which are responsible for the difference between $\delta$ and $I$. Within this approximation, the saturation density $\rho_{sat}$ of asymmetric matter is then still given by equation (6), but replacing $\delta$ by $I$. This no-skin approximation therefore modifies the bulk energy equation (10), and the isoscalar surface energy $E_s^{IS}$. The choice of $I$ instead of $\delta$ to compute the saturation density only slightly worsens the predictive power of the total ETF energy with respect to HF calculations, but the relative weight between bulk and surface energies is modified drastically. In particular, this change of variable switches the sign of the symmetry surface energy [17]. The obvious advantage is that analytical results can be obtained without further approximations than the ones developed in section 2.2, as we now describe in detail.

Replacing $\rho_3(r)$ by $I\rho(r)$ and $\rho_{sat,3}(\delta)$ by $I\rho_{sat}(I)$ in equation (15), allows to express the energy density as a function of $\rho(r)$ only. Thus we can follow the same procedure as in Part I, and analytically integrate the energy density. Making a quadratic expansion in $I$ for the kinetic densities $\tau_3$ gives the following expressions:

$$ E_s^{IV} = C_{\text{surf}}^{IV}(\rho_{\text{sat}}(I), \mathbf{X}_{\text{sky}}^{IV}) \frac{a(A, I)}{r_{\text{sat}}(I)} A^{2/3} I^2 + C_{\text{curv}}^{IV}(\rho_{\text{sat}}(I), \mathbf{X}_{\text{sky}}^{IV}) \left(\frac{a(A, I)}{r_{\text{sat}}(I)}\right)^2 A^{1/3} I^2 \\ + C_{\text{ind}}^{IV}(\rho_{\text{sat}}(I), \mathbf{X}_{\text{sky}}^{IV}) \left(\frac{a(A, I)}{r_{\text{sat}}(I)}\right)^3 I^2 + o\left(\left(\frac{a(A, I)}{r_{\text{sat}}(I)}\right)^4 A^{-1/3} I^2\right), \qquad (23) $$

where $\mathbf{X}_{\text{sky}}^{IV} = \{C_{\text{eff}}, \alpha, D_{\text{eff}}\}$ stands for the effective interaction parameters appearing in the isovector local terms. The coefficients $C_i^{IV}$ are given in appendix B.

As for the isoscalar energy, equation (23) shows that the dominant finite-size effect is a surface term $(\propto A^{2/3})$. Additional finite-size terms, which would be absent in a slab configuration, are found in spherical nuclei. As we have only considered the local part of the isovector energy, we recover the same diffuseness dependence as in the local isoscalar terms equations (17)–(19).

We have seen in section 2.2 that the diffuseness $a$ can be obtained by minimization of the energy per nucleon with respect to its free parameters.

If we neglect the curvature and mass independent terms, we obtain an expression similar to equation (21):
<figure><img src="image_7.png" /><figcaption>Figure 3. As a function of the isospin asymmetry are shown the diffuseness (left side) and total surface energy per nucleon for four isobaric chains (right side, A = 400: full lines, A = 200: dotted lines, A = 100: dashed-dotted lines, A = 50: dashed lines). Red lines: calculations using the slab diffuseness $a_{\text{slab}}$ from equation (24). Blue lines: calculations using the spherical diffuseness $a_{\text{sph}}$ from equation (25). Green line: calculations using the quadratic diffuseness $a_{\text{HF}}$ fitted from HF density profiles [8]. Grey line: calculations using the no-skin diffuseness $a_{\text{TK}}$ from equation (27), based on [20].</figcaption></figure>

$$a_{\text{slab}} = \sqrt{\frac{C_s^{\text{NL}}(I)}{C_{\text{surf}}^{\text{L}}(I) + C_{\text{surf}}^{\text{IV}}(I)I^2}}, \quad (24)$$

where the coefficients $C_{\text{surf}}^i$ depend on the saturation density $\rho_{\text{sat}}(I)$. This expression corresponds to the diffuseness of one-dimensional semi-infinite asymmetric matter. Considering all the terms of equation (23), the diffuseness corresponding to the complete variational problem is given by the solution of the following equation:

$$3(C_{\text{ind}}^{\text{L}} + C_{\text{ind}}^{\text{IV}}I^2)\left(\frac{a_{\text{sph}}}{r_{\text{sat}}}\right)^4 + 2(C_{\text{curv}}^{\text{L}} + C_{\text{curv}}^{\text{IV}}I^2)A^{1/3}\left(\frac{a_{\text{sph}}}{r_{\text{sat}}}\right)^3 \\ + \left(C_{\text{surf}}^{\text{L}} + C_{\text{curv}}^{\text{IV}}I^2\right)A^{2/3} + \frac{1}{r_{\text{sat}}^2}C_{\text{ind}}^{\text{NL}}\left(\frac{a_{\text{sph}}}{r_{\text{sat}}}\right)^2 - \frac{1}{r_{\text{sat}}^2}C_{\text{surf}}^{\text{NL}}A^{2/3} = 0. \quad (25)$$

Figure 3 displays the results of equations (24) and (25). At variance with figure 2 where we only took into account the isoscalar energy, we can clearly see that adding the isovector energy to the variational procedure leads to the expected behavior of a diffuseness increasing with asymmetry. This behavior shows the importance of the isovector part to correctly determine the diffuseness $a$. As for symmetric nuclei, we observe again that the mass dependence of the diffuseness calculated in the spherical case is negligible (only a slight spread of the blue curves, no spread in the red curves). The analytical total surface energy $E_s = E_s^{\text{IS}} + E_s^{\text{IV}}$ per nucleon, given by equations (16) and (23), is plotted on the right side of the same figure.

The results using the slab diffuseness (full red curves) are very close to the ones obtained by solving equation (25) (dash-dotted blue curves), and to the ones using the numerical fit to HF calculations of [8] (dashed green curves), even if the corresponding values for the $a$ parameter are very different. We can see that curvature (and mass-independent) terms are not crucial to determine the diffuseness, though they are important to reproduce the energetics. Therefore the diffuseness can be well determined by the simplest expression, equation (24).
For completeness, we also compare our results to the approximation for the surface energy proposed in [20], and represented by grey curves in figure 3:

$$E_s = E_s^{\text{IS}}(I=0) + 2 \left[ \frac{E_s^{\text{IS}}(I=0)}{A^{2/3}} \frac{L_{\text{sym}}}{K_{\text{sat}}} - \frac{a \left( L_{\text{sym}} - \frac{K_{\text{sym}}}{12} \right)}{r_{\text{sat}}(I=0)} \right] A^{2/3} I^2. \quad (26)$$

In [20], no expression for the diffuseness was proposed. For consistency, we have determined the a parameter entering equation (26) by minimizing the surface energy given by the same equation, leading to:

$$a_{TK} = \sqrt{\frac{C_{\text{surf}}^{\text{NL}}(I=0) \left(1 + \frac{2L_{\text{sym}}}{K_{\text{sat}}} I^2\right)}{C_{\text{surf}}^{\text{L}}(I=0) \left(1 + \frac{2L_{\text{sym}}}{K_{\text{sat}}} I^2\right) - 2 \left(L_{\text{sym}} - \frac{K_{\text{sym}}}{12}\right) I^2}}. \quad (27)$$

To obtain equation (26), the authors of [20] did the same approximation $\rho_3(r) = I\rho(r)$ as we made, neglected the curvature and constant terms, and assumed the equality $E_s^L = E_s^{NL}$ for the isovector part in order to evaluate the non-local isovector energy. As it is shown in Part I, this property fails in a three-dimensional system. As a consequence, the diffuseness which is determined by the balance between local and non-local parts, is overestimated and finally leads to a largely underestimated energy, as seen in figure 3.

To check the accuracy of our analytical no-skin expression given by equations (16), (23) and (24), we will quantitatively compare our analytical results with HF calculations in section 3.3.

## 3.2. Gaussian approximation

To take into account isospin inhomogeneities, we develop in this section an alternative Gaussian approximation of the isovector surface energy. In particular, as in section 2, we will distinguish the bulk asymmetry $\delta$ equation (9) from the global one, $I$, which allows the consideration of skin and Coulomb effects. This approximation is therefore expected to be more realistic than the no-skin procedure developed in section 3.1.

Since $E_s^{IV}$ is the surface isovector energy, the corresponding energy density

$$\mathcal{H}_s^{\text{IV}}[\rho, \rho_3] = \mathcal{H}^{\text{IV}}[\rho, \rho_3] - \frac{\mathcal{H}^{\text{IV}}[\rho_{\text{sat}}, \rho_{\text{sat}3}]}{\rho_{\text{sat}}(\delta)} \rho \quad (28)$$

is negligible at the nucleus center, where $\rho \to \rho_{\text{sat}}$. This is shown in figure 4, which displays this quantity for several nuclei in a representative calculation using the diffusenesses $a$ and $a_p$ from [8], and with the interaction SLy4. Moreover, as it is a surface energy, the maximum is expected to be close to the surface radius $R$, which is the inflection point where $\rho(R) = \rho_{\text{sat}}(\delta)/2$. Thus we approximate the isovector energy density by a Gaussian peak at $r = R$:

$$\mathcal{H}_s^{\text{IV}}(r) \approx \mathcal{G}_{\text{tot}}(r) = \mathcal{A}(A, \delta) e^{-\frac{(r-R)^2}{2\sigma^2(A, \delta)}}, \quad (29)$$

where $\mathcal{A}$ is the maximum amplitude of the Gaussian and $\sigma^2$ its variance at $R$:

$$\mathcal{A}(A, \delta) = \mathcal{H}_s^{\text{IV}}[\rho(R), \rho_3(R)]; \quad \sigma^2(A, \delta) = -\mathcal{A}(A, \delta) \left( \frac{d^2 \mathcal{H}_s^{\text{IV}}}{dr^2} \right)_{r=R}^{-1}. \quad (30)$$
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Numerical isovector energy density profile (red full lines) and Gaussian approximation equation (29) (black dashed-dotted lines) for two masses A = 50 (left curves of each panel) and A = 200 (right curves of each panel). (a) $\delta = 0.1$; (b) $\delta = 0.2$; (c) $\delta = 0.3$; (d) $\delta = 0.4$.</figcaption></figure>

Figure 4 shows the quality of this Gaussian approximation on the energy density profile for several nuclei. Each panel corresponds to a different representative value of $\delta$: $\delta = 0.1$ (upper left) corresponds to most stable nuclei (see figure 1); medium-heavy neutron rich nuclei synthesized in modern radioactive ion facilities lay around $\delta = 0.2$ (upper right); the (largely unexplored) neutron drip-line closely corresponds to $\delta = 0.3$ (lower left); the higher value $\delta = 0.4$ (lower right) is only obtained beyond the dripline, that is for nuclei which are in equilibrium with a neutron gas in the inner crust of neutron stars. We can see that for all these very different asymmetries, the exact energy density (full red lines) is indeed peaked at the equivalent hard sphere radius R. However, we can notice that the profiles have small negative components. We thus expect the Gaussian approximation will overestimate the isovector energy part.

As Gaussian functions and their moments are analytically integrable, this approximation allows obtaining an analytical expression for the isovector energy $E_s^{IV} \approx 4\pi \int r^2 G_{tot}(r) dr$. Indeed, neglecting the terms $\sim e^{-R^2/(2\sigma^2)}$, we obtain (see appendix A):

$$E_s^{IV} = 2(2\pi)^{3/2}\sigma A r_{sat}^2(\delta) \left[ A^{2/3} + \frac{\sigma^2}{r_{sat}^2(\delta)} - \frac{2\pi^2}{3} \left( \frac{a(A, \delta)}{r_{sat}(\delta)} \right)^2 \right], \quad (31)$$

where both $\sigma$ and $A$ depend on the nuclear mass number A, bulk asymmetry $\delta$, and diffusenesses $a$ and $a_p$. The neglected terms are of the order $(a/r_{sat})^4 A^{-2/3}$. We can notice that the curvature term $(\propto A^{1/3})$ is missing. This is due to our approximation. Indeed, we have assumed that the isovector energy is symmetric with respect to the inflection point for which the curvature is zero, such that the curvature is disregarded by construction.

Though equation (31) is an analytical expression, the explicit derivation of the amplitude $A(A, \delta)$ and of the variance $\sigma(A, \delta)$ leads to formulae that are far from being transparent. In particular, it is not clear how the different physical ingredients of the energy functional (compressibility, effective mass, symmetry energy) and of the nucleus properties (neutron skin, diffuseness) affect the isovector surface properties. For this reason, we turn to develop a further approximation for the isovector energy part $E_s^{IV}$ in terms of the nuclear matter
coefficients J, L and K, and of the neutron skin thickness. Moreover, these approximations will allow to find a simple analytical expression for the diffuseness. Making the usual quadratic assumption for the symmetry energy $H^{IV}[\rho, \rho_3] = H_{\text{sym}}[\rho](\rho_3/\rho)^2$, the amplitude $A(A, \delta)$ equation (30) reads

$$A(A, \delta) = H_{\text{sym}}[\rho(R)] \left( \frac{\rho_3(R)}{\rho(R)} \right)^2 - H_{\text{sym}}[\rho_{\text{sat}}(\delta)] \frac{\rho(R)}{\rho_{\text{sat}}(\delta)} \delta^2. \quad (32)$$

In order to have a simpler explicit expression, we make a density expansion of the symmetry energy per nucleon $e_{\text{sym}}[\rho] = H_{\text{sym}}[\rho]/\rho$ around a density $\rho_*$, such that:

$$H_{\text{sym}}[\rho] = \rho \left[ J_* + \frac{L_*}{3\rho_*}(\rho - \rho_*) + \frac{K_*}{18\rho_*^2}(\rho - \rho_*)^2 \right], \quad (33)$$

where $J_* = H_{\text{sym}}[\rho_*/\rho_*, L_* = 3\rho_*\partial(H_{\text{sym}}/\rho_*)/\partial\rho|_{\rho_*}$ and $K_* = 9\rho_*^2\partial^2(H_{\text{sym}}/\rho_*)/\partial\rho^2|_{\rho_*}$. As we can see in equation (32), we need to evaluate the symmetry energy at two different densities: at $\rho_{\text{sat}}(\delta)$ and at the surface radius where $\rho(R) = \rho_{\text{sat}}(\delta)/2$. For this reason, we will apply equation (33) to two different densities $\rho_* = \rho_{\text{sat}}(0)$ and $\rho_* = \rho_{\text{sat}}(0)/2$. At $\rho_* = \rho_{\text{sat}}(0)$, the coefficients $(J_*, L_*, K_*)$ are the usual symmetry energy coefficients $(J_{\text{sym}}, L_{\text{sym}}, K_{\text{sym}})$. At one half of the saturation of symmetric nuclear matter, $\rho_* = \rho_{\text{sat}}(0)/2$ we label the corresponding coefficients $(J_{1/2}, L_{1/2}, K_{1/2})$. Using the expansion around $\rho_* = \rho_{\text{sat}}(0)/2$ for the first term of equation (32) and around $\rho_* = \rho_{\text{sat}}(0)$ for the second one, we obtain, at second order in $\delta$:

$$\frac{\mathcal{A}(A, \delta)}{\rho_{\text{sat}}(0)} = \frac{J_{1/2}}{8} \left(\frac{\Delta R(a)}{a(A, \delta)}\right)^2 + \frac{J_{1/2}}{2} \left[ \frac{\Delta R(a)}{a(A, \delta)} - \frac{1}{2} \left(\frac{\Delta R(a)}{a(A, \delta)}\right)^2 \right] \delta + \frac{J_{1/2}}{2} \left[ \left(1 - \frac{J_{\text{sym}}}{J_{1/2}}\right) - \frac{\Delta R(a)}{a(A, \delta)} \frac{1}{4} \left(1 + \frac{L_{\text{sym}}L_{1/2}}{J_{1/2}K_{\text{sat}}}\right) \left(\frac{\Delta R(a)}{a(A, \delta)}\right)^2 \right] \delta^2. \quad (34)$$

Notice that the $K_{\text{sym}}$ parameter does not appear in this equation because of the truncation at second order in $\delta$. In equation (34), the isospin asymmetry inhomogeneities clearly appear through the quantity $\Delta R(a) = R(a) - R_p(a)$ which represents the neutron skin thickness:

$$\Delta R(a) = \Delta R_{\text{HS}} \left( 1 + \frac{\pi^2}{3} \frac{a^2}{R_{\text{HS}} R_{\text{HS},p}} \right), \quad (35)$$

where $\Delta R_{\text{HS}}(A, Z) = \Delta R(a = 0, A, Z) = R_{\text{HS}}(A) - R_{\text{HS},p}(Z)$ is the neutron skin thickness of nuclei theoretically described by hard spheres. Moreover, we have considered the diffuseness difference $a - a_p$ as a second-order correction with respect to the neutron skin, and have assumed $a = a_p$ in equation (34). We have also used the following expansion in $\Delta R(a)/a$ to evaluate $\rho_3(R)$:

$$2\rho_p(R) = \rho_{\text{sat},p}(\delta)[1 - \Delta R(a)/(2a)] + O((\Delta R(a)/a)^3).$$

Equation(34) gives a relatively simple and transparent expression of the isovector energy density at the nuclear surface, as a function of the equation of state parameters. The situation is more complicated for the variance $\sigma(A, \delta)$ which also enters the isovector energy equation (31). This quantity involves the second spatial derivative of the energy density equation (30), therefore its explicit expression is not transparent, even with the previous simplifications. Extra approximations are in order.
From figure 4, we can observe that the width of the numerical Gaussians, that is the values of $\sigma^2(A, \delta)$, is almost independent of the bulk isospin $\delta$. This numerical evidence can be understood from the fact that the width gives a measure of the nucleus surface, which is mostly determined by isoscalar properties. It is therefore not surprising that the dominant isospin dependence is given by the amplitude $A$ which represents the isovector energy density at the surface. For this reason, we evaluate the variance at $\delta = 0$:

$$ \sigma(A, \delta) \approx \sigma(A) = \sqrt{\frac{2}{1 - \frac{K_{1/2}}{18J_{1/2}}}} a_0 = \sigma_0. \qquad (36) $$

In this equation, $a_0$ stands for the diffuseness at $\delta = 0$. We recall that this quantity does not depend on the nucleus mass if we do not take into account terms beyond surface in the variational derivation. This approximate mass independence of the variance can be verified in figure 4: the width of the two Gaussians corresponding to $A = 50$ and $A = 200$ are very close. Neglecting the isovector component at $\delta = 0$, the diffuseness is then given by the expression (21) valid for symmetric matter:

$$ a_0 = \sqrt{C_{\text{surf}}^{\text{NL}}(\delta = 0)/C_{\text{surf}}^{\text{L}}(\delta = 0)}. \qquad (37) $$

Inserting equations (34) and (36) into (31), the surface isovector energy can be expressed as a function of the symmetry energy coefficients ($J_{\text{sym}}$, $L_{\text{sym}}$, $K_{\text{sym}}$):

$$ \begin{aligned} E_s^{\text{IV}} = {}& 3 \sqrt{\frac{\pi}{1 - \frac{K_{1/2}}{18J_{1/2}}}} \frac{\rho_{\text{sat}}(0)}{\rho_{\text{sat}}(\delta)} \frac{a_0}{r_{\text{sat}}(\delta)} J_{1/2} \\ & \times \left\{ \frac{1}{4} \left( \frac{\Delta R(a)}{a(A, \delta)} \right)^2 + \left[ \frac{\Delta R(a)}{a(A, \delta)} - \frac{1}{2} \left( \frac{\Delta R(a)}{a(A, \delta)} \right)^2 \right] \delta \right. \\ & \left. + \left[ \left( 1 - \frac{J_{\text{sym}}}{J_{1/2}} \right) - \frac{\Delta R(a)}{a(A, \delta)} - \frac{1}{4} \left( 1 + \frac{L_{\text{sym}} L_{1/2}}{J_{1/2} K_{\text{sat}}} \right) \left( \frac{\Delta R(a)}{a(A, \delta)} \right)^2 \right] \delta^2 \right\} \\ & \times \left\{ A^{2/3} + \frac{2}{1 - \frac{K_{1/2}}{18J_{1/2}}} \left( \frac{a_0}{r_{\text{sat}}(\delta)} \right)^2 - \frac{2\pi^2}{3} \left( \frac{a(A, \delta)}{r_{\text{sat}}(\delta)} \right)^2 \right\}. \end{aligned} \qquad (38) $$

In principle the surface coefficients ($J_{1/2}$, $L_{1/2}$, $K_{1/2}$) can be expressed as a function of the bulk ones ($J_{\text{sym}}$, $L_{\text{sym}}$, $K_{\text{sym}}$) by using polynomial expansion in the density. However, we can see from equation (38) that the surface isovector energy $E_s^{\text{IV}}$ is proportional to the symmetry energy $J_{1/2}$ evaluated at the surface $R$. It is quite natural that the surface energy component is mainly determined by the surface properties of the nuclei, and therefore, the surface symmetry energy is mainly proportional to the isovector parameter $J_{1/2}$. For this reason, expressing equation (38) only in terms of bulk quantities ($J_{\text{sym}}$, $L_{\text{sym}}$, $K_{\text{sym}}$) would make equation (38) less transparent.

For completely symmetric nuclei, that is $\Delta R = 0$ and $\delta = 0$, the isovector energy is identically zero as it should. However, if we neglect the neutron skin thickness only, that is we consider $\Delta R = 0$ but $\delta \neq 0$, a non-zero isovector surface energy is obtained, given by
$$E_{\text{surf}}^{\text{IV}, \Delta R=0} = 3 \sqrt{\frac{\pi}{1 - \frac{K_{1/2}}{18J_{1/2}}}} \frac{\rho_{\text{sat}}(0)}{\rho_{\text{sat}}(\delta)} \frac{a_0}{r_{\text{sat}}(\delta)} (J_{1/2} - J_{\text{sym}}) \delta^2 A^{2/3}. \quad (39)$$

This expression is proportional to the energy density difference between bulk and surface ($J_{1/2} - J_{\text{sym}}$), that is to the $L_{\text{sym}}$ parameter. In this approximation, the diffuseness $a(A, \delta)$ does not appear, which means that the isovector surface energy contributes to the determination of the diffuseness only if we consider the neutron skin.

From a mathematical point of view we can also consider the limit $\delta = 0$, $\Delta R \neq 0$, giving:

$$E_{\text{surf}}^{\text{IV}, \delta=0} = \frac{3}{4} \sqrt{\frac{\pi}{1 - \frac{K_{1/2}}{18J_{1/2}}}} J_{1/2} \frac{\Delta R^2(a_0)}{a_0 r_{\text{sat}}(0)} A^{2/3}. \quad (40)$$

This expression shows that an isovector surface energy can be induced in asymmetric nuclei even if no asymmetry is present in the bulk. Of course in realistic situations the bulk asymmetry and the difference between neutron and proton radii are not independent variables; in particular the skin is negligible if $\delta = 0$ as we have already assumed in order to obtain equation (37) above.

Equation (38) shows that even in our rather crude approximation the surface symmetry energy presents a very complex dependence on the physical quantities that measure isospin inhomogeneity, namely the bulk asymmetry $\delta$ and the neutron skin thickness $\Delta R$. In particular we find that $E_s^{\text{IV}}(A, \delta)$ is not quadratic with $\delta$ but has non-negligible linear components (see also figure 7 below). We have also quantitatively tested that both linear and quadratic terms in $\Delta R$ are required to correctly reproduce the surface isovector energy. It is interesting to notice that the linear components mix $\delta$ and $\Delta R$. Indeed, as we can see in equations (39) and (40), putting to zero one of those variables, which both measure the isospin inhomogeneities, leads to a quadratic behavior with respect to the other variable (cf. equations (39) and (40)).

Similar to the previous section, the diffuseness is the only unconstrained parameter of the model. It can therefore be determined in a variational approach by minimizing the total (isoscalar and isovector) surface energy. We have already discussed the fact that only the dominant $\propto A^{2/3}$ terms are important to evaluate the diffuseness. For this reason, we neglect again terms beyond plane surface, and we approximate the neutron skin thickness $\Delta R$ by the hard sphere approximation $\Delta R_{\text{HS}}$. Neglecting the quadratic terms in the expansion in $\Delta R_{\text{HS}}/a$, we obtain

$$a_{\text{Gauss}}^2(A, \delta) = a_{\text{IS}}^2(\delta) + \sqrt{\frac{\pi}{1 - \frac{K_{1/2}}{18J_{1/2}}}} \frac{\rho_{\text{sat}}(0)}{\rho_{\text{sat}}(\delta)} \frac{3J_{1/2}(\delta - \delta^2)}{C_{\text{surf}}^{\text{L}}(\delta)} a_0 \Delta R_{\text{HS}}(A, \delta), \quad (41)$$

where $a_{\text{IS}}(\delta)$ is the diffuseness obtained by neglecting the isovector component: $a_{\text{IS}}(\delta) = \sqrt{C_{\text{surf}}^{\text{NL}}(\delta)/C_{\text{surf}}^{\text{L}}(\delta)}$. We found in section 2 that $a_{\text{IS}}$ slightly decreases with the isospin asymmetry (see figure 2), which does not appear to be consistent with the behavior observed in full HF calculations. Now considering in the variational principle the isovector term in addition to the isoscalar one, the diffuseness $a$ given by equation (41) acquires an additional term which modifies its global $\delta$ dependence. The complete result equation (41) is displayed in figure 5. We can see that the additional term due to the isovector energy contribution inverses the trend found section 2, as expected. More specifically, though it does not clearly appear in equation (41), the analytical diffuseness is seen to quadratically increase with $\delta$, corroborating the assumption found in [8].
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Diffuseness (left side) and surface energy per nucleon for four different isobaric chains (right side), as a function of the bulk isospin asymmetry. Red lines: $a_{\text{gaus}}$ from equation (41). Blue lines: $a^{\text{min}}$, $a_p^{\text{min}}$ from the minimization of the exact numerically calculated ETF surface energy. Green lines: total $a_{\text{HF}}$ and proton $a_{p,\text{HF}}$ diffuseness fitted from HF density profiles, taken from [8].</figcaption></figure>

Although we only considered terms $\propto A^{2/3}$, as in a slab geometry, the results slightly depend on the nucleus mass as shown by the slight dispersion of the different red curves in figure 5. This is due to the neutron skin since $\Delta R_{\text{HS}}(A, \delta)$ increases with decreasing mass number A. For comparison, the diffusenesses $a_{\text{HF}}$ and $a_{p,\text{HF}} \neq a$ obtained by a fit of HF density profiles in [8] are also represented in figure 5 (green curves), as well as the numerically calculated pair ($a^{\text{min}}$, $a_p^{\text{min}}$) which minimises the energy (blue curves).

As we can see, these diffusenesses differ from each other significantly, but their consequence on the energy is small as we can observe in the right part of figure 5, which displays the corresponding surface energy $E_s = E_s^{\text{IS}} + E_s^{\text{IV}}$ per nucleon, for different isobaric chains. In this figure, the blue curves correspond to a numerical integration of the ETF energy density, using the diffusenesses which minimize the total surface energy. These results can thus be considered as 'exact' ETF results. The use of the very different $a$ and $a_p$ values fitted from HF (green lines) leads to only slightly different energies, except for the lightest isobar chain. The analytical approximation given by the sum of equations (16) and (38), is also plotted (red curves), where the diffuseness is given by the analytical formula equation (41). We can see that our analytical approximation closely follows the 'exact' ETF results. All the curves show a positive surface symmetry energy, which contrasts with figure 3. As it has been discussed in [17], this change of sign is due to the choice between the bulk asymmetry $\delta$ or the global asymmetry $I$, in the definition of the bulk energy. This choice obviously affects the residual part of the energy $E_s$, since the sum of the two gives the same ETF functional. This residual part is, to first order, given by the surface symmetry energy as discussed in [17].

In order to further validate the analytical results of this section, quantitative comparisons with HF calculations are shown in the next section 3.3.

### 3.3. Comparison to HF calculations

In this section, we explore the level of accuracy of both the no-skin approximation and the Gaussian approximation, respectively, developed in sections 3.1 and 3.2.

As previously discussed, the two different approximations lead to two different bulk energetics. Neglecting isospin inhomogeneities implies that the bulk asymmetry $\delta$ is equalized to the average asymmetry $I$. Thus the bulk quantities $\rho_{\text{sat}}$ and $E_b$ defined by equations (6) and
<figure><img src="image_10.png" /><figcaption>Figure 6. Total energy $E = E_b + E_s$ per nucleon (upper part) and difference between the nuclear HF energy per nucleon and the prediction of the different models (lower part) as a function of the nucleus asymmetry $I = 1 - 2Z/A$ calculated within the no-skin approximation, equation (42) (blue dotted lines) within the Gaussian approximation, equation (43) (red full lines), and from the model in [20] (green dashed lines). Nuclear HF energy is given by the stars. (a) $A = 50$; (b) $A = 100$; (c) $A = 200$; (d) $A = 400$.</figcaption></figure>

(10) depend on $I$, and the total energy of a nucleus ($A, I$) within the no-skin approximation is given by

$$E_{\text{NoSkin}}(A, I) = E_b(A, I) + E_s^{\text{IS}}(A, I) + E_s^{\text{IV}}(A, I), \quad (42)$$

where $E_s^{\text{IS}}(A, I)$ is given by equation (16) (with $I$ instead of $\delta$), $E_s^{\text{IV}}(A, I)$ by equation (23), and the diffuseness is given by equation (24).

On the other hand, the Gaussian approximation allows defining two independent density profiles. Therefore, the bulk energy depends on the bulk asymmetry $\delta(A, I)$ defined by equation (9) and the total energy of a nucleus ($A, I$) within this approximation is given by
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Bulk (upper left), surface (lower left), isoscalar (upper right) and isovector (lower right) surface energy per nucleon as a function of the bulk asymmetry $\delta$ for isobaric nuclei $A = 100$, predicted by equation (43). Different Skyrme interactions are considered: SLy4 [27] (full red), SkI3 [37] (dashed green), SGI [38] (dotted blue), LNS [39] (dashed-dotted black).</figcaption></figure>

$$E_{\text{Gauss}}(A, I) = E_b(A, \delta) + E_s^{\text{IS}}(A, \delta) + E_s^{\text{IV}}(A, \delta), \quad (43)$$

where $E_s^{\text{IS}}(A, \delta)$ is given by equation (16), $E_s^{\text{IV}}(A, \delta)$ by equation (38), and the isoscalar diffuseness is given by equation (41).

In figure 6, we compare the analytical expressions (42) and (43) with HF energy calculations for different isobaric chains. To compare the same quantities, we used the same interaction (SLy4) and removed the Coulomb energy from the HF calculation. We can see from the figure that the no-skin and the Gaussian approximations predict close values for the total energy. For low asymmetries $I \lesssim 0.2$ where the two models are almost indistinguishable: they reproduce the microscopic calculations with a very good accuracy, especially for medium-heavy nuclei $A \gtrsim 100$. However, for higher asymmetries $I \gtrsim 0.2$ where the symmetry energy becomes important, a systematic difference between the two models appears and increases up to $\sim 400 \text{ keV A}^{-1}$ for the highest asymmetries $I \sim 0.4$: the Gaussian approximation is systematically closer to the microscopic results than the no-skin model. This observation highlights the importance of taking into account the isospin asymmetry inhomogeneities, considering the neutron skin and at the same time differentiating the bulk asymmetry $\delta$ from the global one $I$, as was discussed in [17]. For comparison, the predictions of [20] are also represented (dashed green line and square symbols). We recall that in that work the isospin inhomogeneities were neglected in a similar way to our no-skin approximation. In addition to that, the approximation was made that the non-local part of the surface energy is equal to the local part. This approximation is already not justified in $N = Z$ nuclei, as was discussed at length in Part I. Moreover, it turns out that it becomes increasingly false with increasing isospin asymmetry (see figure 9, below), and leads to an isospin dependence for the diffusivity parameter which is not compatible with the results of HF calculations (see
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Diffuseness $a_{\text{Gauss}}$ (upper panel), neutron skin thickness $\Delta R$ (middle panel) and the ratio $\Delta R/a_{\text{Gauss}}$ (lower panel) as a function of the bulk asymmetry $\delta$ for the isobaric chain A = 100, predicted within the Gaussian approximation (see text). Different Skyrme interactions are considered: SLy4 [27] (full red), SkI3 [37] (dashed green), SGI [38] (dotted blue), LNS [39] (dashed-dotted black).</figcaption></figure>

figure 3, above). As a consequence, the corresponding analytical mass formula strongly deviates from the microscopic calculation for $I > 0.2$. Concerning the formulas proposed in the present work, the accuracy of equation (43) is better than ~200 keV $A^{-1}$ for medium-heavy nuclei, which is similar to the predictive power of the bulk part of spherical HF calculations for this effective interaction, with respect to experimental data.

To conclude, the Gaussian approximation developed in section 3.2 provides a reliable analytical formula, especially for the surface symmetry energy. For this reason we will only use the Gaussian approximation to further study the different components of the nuclei energetics, as we turn to do in the next section.

**4. Study of the different energy terms**

In this section, we use the analytical formulae based on the Gaussian approximation detailed in section 3.2, to study the different components of nuclear energetics. As we have previously discussed throughout this paper, we can decompose the nucleus total energy $E$ into bulk $E_b$ and surface $E_s$ parts. Both can be written as sums of isoscalar $E_i^{\text{IS}}$, that is the part independent of $\rho_3(r)$, and isovector $E_i^{\text{IV}}$ terms. The surface energy can be further split into plane surface
<figure><img src="image_13.png" /><figcaption>Figure 9. Decomposition of the local (left) and non-local (right) part of the isoscalar surface energy per nucleon, into its surface and curvature component as a function of the bulk asymmetry δ for the isobaric chain A = 100, as predicted by equation (43). Different Skyrme interactions are considered: SLy4 [27] (full red), SkI3 [37] (dashed green), SGI [38] (dotted blue), LNS [39] (dashed-dotted black).</figcaption></figure>

$E_{surf} \propto A^{2/3}$, curvature $E_{curv} \propto A^{1/3}$ and mass independent $E_{ind}$ terms. Finally, we can distinguish the local $E_{i}^{\text{IS,L}}$ and the non-local $E_{i}^{\text{IS,NL}}$ components of the surface isoscalar part only, since we did not discriminate them in the Gaussian approximation used for the isovector energy. In summary, the energy of an (A, I) nucleus can be written as

$$E(A, I) = E_b(A, \delta) + E_s(A, \delta), \qquad (44)$$

$$E_s(A, \delta) = E_s^{\text{IS}}(A, \delta) + E_s^{\text{IV}}(A, \delta), \qquad (45)$$

$$E_s^{\text{IV}}(A, \delta) = E_{\text{surf}}^{\text{IV}}(A, \delta) + E_{\text{ind}}^{\text{IV}}(A, \delta), \qquad (46)$$

$$E_s^{\text{IS}}(A, \delta) = E_{\text{surf}}^{\text{IS}}(A, \delta) + E_{\text{curv}}^{\text{IS}}(A, \delta) + E_{\text{ind}}^{\text{IS}}(A, \delta), \qquad (47)$$

$$E_{\text{surf}}^{\text{IS}}(A, \delta) = E_{\text{surf}}^{\text{IS,L}}(A, \delta) + E_{\text{surf}}^{\text{IS,NL}}(A, \delta), \qquad (48)$$

$$E_{\text{curv}}^{\text{IS}}(A, \delta) = E_{\text{curv}}^{\text{IS,L}}(A, \delta) + E_{\text{curv}}^{\text{IS,NL}}(A, \delta), \qquad (49)$$

where the bijective relation (for a given mass) between I and δ is given by equation (9). The different isoscalar terms $E_{i}^{\text{IS,j}}$ are defined by equations (16) to (19), with the diffuseness $a(A, \delta)$ determined within the Gaussian approximation, equation (41). The isovector components $E_{i}^{\text{IV}}$ are introduced in equation (38), where the curvature term, in this Gaussian approximation, is identically zero by construction.

In the following, we will study each of these terms, and specifically their dependence with the asymmetry δ. For this comparison, we have chosen a representative isobaric chain A = 100 for which the ETF approximation was successfully compared to HF results in figure 6, for the SLy4 interaction. For this choice of mass, $δ \approx 0$ corresponds to the proton dripline and $δ \approx 0.3$ the neutron dripline (see figure 1).

Due to our limited experimental knowledge of the isovector properties of the effective interaction, the behavior of the different energy terms with asymmetry is to some extent model-dependent. In order to sort out general trends we have considered different Skyrme functionals which approximately span the current uncertainties on the density dependence of
<table><caption>**Table 1.** Bulk and surface nuclear properties for the different Skyrme interactions.</caption><thead><tr><th>Interaction</th><th>ρ<sub>sat</sub> (0)<br>(fm<sup>-3</sup>)</th><th>m*/m</th><th>J<sub>sym</sub><br>(MeV)</th><th>L<sub>sym</sub><br>(MeV)</th><th>K<sub>sym</sub><br>(MeV)</th><th>J<sub>1/2</sub><br>(MeV)</th><th>L<sub>1/2</sub><br>(MeV)</th><th>K<sub>1/2</sub><br>(MeV)</th></tr></thead><tbody><tr><td>SLY4 [27]</td><td>0.1595</td><td>0.595</td><td>32.00</td><td>46.0</td><td>-119.8</td><td>22.13</td><td>38.6</td><td>-74.0</td></tr><tr><td>SkI3 [37]</td><td>0.1577</td><td>0.577</td><td>34.83</td><td>100.5</td><td>73.0</td><td>18.85</td><td>46.7</td><td>-25.2</td></tr><tr><td>SGI [38]</td><td>0.1544</td><td>0.608</td><td>28.33</td><td>63.9</td><td>-52.0</td><td>16.75</td><td>38.4</td><td>-29.7</td></tr><tr><td>LNS [39]</td><td>0.1746</td><td>0.826</td><td>33.43</td><td>61.5</td><td>-127.4</td><td>21.10</td><td>44.6</td><td>-56.8</td></tr></tbody></table>

the symmetry energy. The corresponding bulk parameters are reported in table 1. In this table, the calculated surface coefficients ($J_{1/2}$, $L_{1/2}$, $K_{1/2}$) entering equations (38) and (41) are also given. As is well known [36], the different interactions are very close at half saturation density, reflecting the fact that all Skyrme parameters have been fitted on ground state properties of finite nuclei, which correspond to an average density of the order of $\rho_{sat}/2$. Nevertheless, a considerable spread is already seen at saturation density, showing that the extrapolation of isovector properties to unexplored density domains is still not well controlled [36]. Concerning the LNS interaction, the parametrization proposed in [39] corresponds to a too high saturation density which is not realistic. This induces a trivial deviation with respect to the other interactions in both the bulk and surface isovector components. For this reason, only the isovector properties of this functional are of interest for this study.

A more complete study of the effective interactions parameter space would be necessary to reach sound conclusions on the quantitative model dependence, but from the representative chosen interactions, we can already address some qualitative interpretations.

The bulk energy per nucleon is shown in the upper left panel of figure 7. At low asymmetries, the curves are indistinguishable, reflecting the good present knowledge of symmetric nuclear matter properties. The only exception is given by LNS, which presents a global shift with respect to the other functionals. As already remarked, this is due to the irrealistically high saturation density of this parametrization (table 1). However, we can see that the behavior with isospin is comparable to the one of the other functionals, reflecting a compatible bulk symmetry energy. For the highest asymmetries $\delta \gtrsim 0.25$, we can see that all the parametrizations differ, which reflects the larger uncertainties for asymmetric matter.

The lower left panel of figure 7 displays the surface corrections. We can see that the qualitative behavior of the different models is the same: $E_s/A$ increases with the asymmetry, leading to a positive sign of the corresponding symmetry energy. As has been already discussed in [17], this comes from the consideration of the bulk asymmetry $\delta$ instead of the global one, $I$, in the definition of the nuclear bulk. The increase rate with isospin is not the same in the different models, reflecting the different surface symmetry energies of the functionals. In particular, the steep behavior predicted by the SkI3 parametrization is due to the stiff isovector properties of this effective interaction (see $L_{sym}$ and $K_{sym}$ in table 1), which lie close to the higher border of the presently accepted values for these parameters [36].

The right part of figure 7 shows the energy decomposition of equations (45) and (46). As expected, at $\delta = 0$, though not identically zero (see equation (40)), the isovector energy (lower panel) is completely negligible. This *a posteriori* justifies the assumption $E_s^{IV}(0) = 0$ we made in order to obtain $a_0$ in equation (36). However, for asymmetric systems, though smaller than the isoscalar energy (upper panel), the isovector energy cannot be neglected. Indeed, its dependence with $\delta$ is much stronger, meaning that the isovector term is the most important term determining the surface symmetry energy. Concerning the mass-independent term, we can see that it is negligible compared to the other components, as expected for the
medium-heavy nucleus concerned by this picture. Finally, we can observe that the isovector energy is not quadratic with $\delta$, thus confirming that the linear terms of equation (38) cannot be neglected.

Concerning the predictions of the different functionals, we can observe that the shift in the bulk properties of LNS with respect to the other models is compensated by an opposite shift in the surface energy. Concerning the three other models, the deviation at high asymmetry is due to the existing uncertainties in the isovector properties of the functionals. As expected, the deviation is the largest in the isovector part (~53% difference between SLy4 and SkI3 at the highest asymmetry) with respect to the isoscalar one (~11% difference). Much less expected is the fact that at $\delta = 0$, where the SLy4, SkI3 and SGI models are in perfect agreement on the bulk energy and the isovector energy is negligible, they however differ from ~500 keV per nucleon on the surface energies. We will come back to this surprising result later in this section.

Figure 8 shows the predictions of the different functionals concerning the parameters associated to the density profiles, namely the diffuseness (upper panel), the neutron skin (middle panel) and their ratio (lower panel). We can see that, for a given asymmetry $\delta$, the spread of the diffuseness values given by equation (41) is very important, reflecting the poor knowledge of this quantity. These large uncertainties can be understood considering that the diffuseness does not seem to affect the energy in a systematic way. In particular, though SkI3 and SGI models surprisingly give the same diffuseness, the corresponding surface properties systematically differ. Moreover, this similarity of the diffuseness cannot be straightforwardly linked to any specific interaction property or parameter (see table 1). This reflects again the fact that the diffuseness is a delicate balance of all energy components, and is determined by very subtle competing and opposite effects.

The middle part of the figure shows the obvious correlation between $\Delta R = R - R_p$ and $\delta$. It is clear from this behavior that quadratic terms in the neutron thickness cannot be neglected to correctly estimate the symmetry energy (see equation (38)). It is interesting to observe that the SGI and LNS models give very close results for this quantity, and the same was true for the isovector part of the surface energy in figure 7 above. This comes from the fact, already observed in the literature [36], that $\Delta R$ is mainly determined by the slope of the symmetry energy $L_{\text{sym}}$ [36], which are close in the SGI and LNS models. Our work confirms that the neutron thickness can be viewed as a measurement of the $L$ parameter. Indeed, $\Delta R$ can be well approximated using the equivalent hard spheres radii $R_{\text{HS}}(\delta)$, $R_{\text{HS},p}(\delta)$, see equation (35). This means that $\Delta R$ can be seen as a function of the saturation density $\rho_{\text{sat}}(\delta)$. In turn, the saturation density is given by equation (6) which at first order is quadratic in $\delta^2$ with the coefficient $L_{\text{sym}}/K_{\text{sat}}$. Since $K_{\text{sat}}$ is relatively well constrained, we then understand why $\Delta R$ is mainly determined by $L_{\text{sym}}$. In particular, the neutron skin thickness is predicted to be the same in the two specific interactions SGI and LNS. Since the surface isovector energy equation (38) at a given bulk asymmetry mainly depends on the neutron skin, this also explains why we obtain the same energies for the two models in figure 7. This essential role of $\Delta R$ to determine the symmetry energy is confirmed observing from figures 7 and 8 that Skyrme models which predict thicker neutron skin, that is higher $L_{\text{sym}}$, give systematically larger values of the isovector surface energy.

The lower part of figure 8 shows the ratio $\Delta R/a$ as a function of $\delta$. Though it is the quantity which mainly governs the behavior of equation (38), it does not constrain the surface isovector energy $E^{\text{IV}}$. Indeed, same $\Delta R/a$ from the functionals SkI3 and SGI lead to different energies (figure 7, lower right panel), corroborating the above discussion: only the $L$ parameter, or equivalently the neutron skin thickness $\Delta R$, is relevant to determine the isovector
contribution. This stresses the importance of the experimental measurement of neutron skin thickness as a key quantity for the knowledge of the density dependence of the symmetry energy [36].

To conclude, we study in figure 9 the decomposition into local and non-local terms as predicted by the different functionals. Only the isoscalar part of the surface energy is considered because these different terms are mixed up in the Gaussian approximation we have employed for the isovector component. Again, we can see that the qualitative behavior of the different Skyrme models is the same for each specific term. We can then safely conclude that the non-local curvature component $E_{curv}^{IS,NL}$ can be neglected for medium-heavy nuclei $A \gtrsim 100$, but the local curvature energy has to be taken into account since it represents for these nuclei 10% to 25% of the total surface local energy, depending on the interaction choice and on the asymmetry $\delta$. Concerning the $\delta$ dependence of the isoscalar surface energies in figure 9, we can notice that the local and non-local parts have opposite behaviors, leading to the rather flat curves observed in figure 7, upper right panel. In section 2.2, we have shown that the exact equality $E_{surf}^{IS,L} = E_{surf}^{IS,NL}$ is obtained only if both curvature and isovector terms are neglected in the determination of the diffuseness (see also Part I). However, the neglect of isovector terms leads to a wrong dependence with $\delta$ as shown in figure 2. Thus, isovector terms cannot be avoided.

The results of figure 9 clearly show that, once these terms are consistently added in the variational procedure (equation (41)), the equality $E_{surf}^{IS,L} = E_{surf}^{IS,NL}$ is completely violated for asymmetric systems. Therefore the isoscalar energy strongly depends on the neutron skin thickness, even if it is an indirect dependence through the diffuseness. This shows that, though the energy can be split into different terms, the latter cannot be decorrelated and they have to be treated together.

We have already observed in figure 7 that the different functionals predict very different surface energy at $\delta = 0$, which might be surprising considering that the symmetric nuclear properties are supposed to be well constrained by experimental data. An obvious interpretation would be that the discrepancy comes from the surface properties, that is the non-local gradient terms and the (poorly constrained) diffuseness parameter. However, comparing the different values of the predicted diffuseness at $\delta = 0$ from figure 8, we can see that $a_{SGI} < a_{LNS} = a_{SkI3} < a_{SLy4}$. This inequality sequence is not respected for the surface energy $E_{surf}(\delta = 0)$ in figure 9, meaning that the difference of surface energies cannot be ascribed to the diffuseness. The possible dependence on the couplings of gradient and spin-orbit terms is also excluded. Indeed, we can see from figure 9 that at $\delta = 0$, the isovector part is zero by definition and therefore the equality $E_{surf}^{IS,L} = E_{surf}^{IS,NL}$ is verified. This means that the total surface energy for symmetric bulk is $E_{surf} = 2E_{surf}^{IS,L}$, which does not depend on the non-local terms of the functional, but only depends on the bulk interaction coefficients $(\rho_{sat}(0), C_0, C_3, C_{eff}, \alpha)$.

We can conclude that the differences of the total surface energies observed for $\delta = 0$, that is nuclei very close to isospin symmetry, in figure 7, does not come from the non-local properties but are intrinsically linked to the bulk interaction coefficients ($C_0, C_3, C_{eff}, \alpha$), though the SLy4, SkI3 and SGI models correspond to compatible isoscalar equations of state (that is: compatible values for the saturation density $\rho_{sat}(0)$, bulk energy $E_b(\delta = 0)$, compressibility $K_{sat}$ and effective mass). This shows that, at variance with the skin thickness $\Delta R$, which is strongly correlated to the isovector equation of state, the nuclear surface energy very poorly constrains the equation of state, even for symmetric or quasi-symmetric nuclei.

The different isoscalar functional parameters ($C_0, C_3, C_{eff}, \alpha$) thus appear largely unconstrained. A detailed correlation study would be very useful to determine the most useful
observables which could be experimentally measured in order to break the observed degeneracy.

**5. Summary and conclusions**

In this paper we have addressed the problem of the determination of an analytical mass formula with coefficients directly linked to the different parameters of standard Skyrme functionals, in the ETF approximation at second order in $\hbar$. The purpose of this effort is twofold. On one side, such a formula is useful for astrophysical applications where extendend calculations are needed covering the whole mass table and using a variety of effective interaction to assess the sensitivity of astrophysical observables to the nuclear physics inputs [21]. On the other side, analytical expressions of the different coefficients of the mass formula in terms of the Skyrme couplings allow a better understanding of the correlation between these couplings and the different aspects of nuclear energetics, for the construction of optimized fitting procedures of the functionals.

The modelling of Fermi density profiles proposed in [8] allows an (almost) exact analytical evaluation of the isoscalar part of the nuclear energy, naturally leading to the appearance in the surface energy of a curvature term and a constant term independent of the baryonic number.

The calculation of the isovector part is highly non-trivial. No exact analytical integration of the ETF functional is possible in the presence of isospin inhomogeneities, and approximations have to be made. We have proposed two different approximations for the determination of the surface symmetry energy. The first approximation consists in completely neglecting the difference between the neutron and proton radius, that is the neutron skin $\Delta R$. The resulting surface energy shows a quadratic dependence on the isospin asymmetry $I$, and consists of local and non-local plane surface, curvature and mass dependent terms which are simple generalizations of the expressions obtained for symmetric nuclei in Part I. Surprisingly, this crude approximation reproduces numerical HF results very well for all stable nuclei up to asymmetries of the order of $I \approx 0.2$, and leads to a relatively limited overestimation of the order of $\approx 400$ KeV/nucleon close to the driplines.

A better approximation is obtained if isospin inhomogeneities are accounted for. To this aim, we have introduced a different radius for the neutron and proton distributions, as well as an explicit difference between the global asymmetry, $I$, and the asymmetry in the nuclear bulk, $\delta$, due to both Coulomb and neutron skin effects. In this more general case, to obtain a mass formula we make the assumption that the surface energy density is peaked at the nuclear surface, and curvature terms can be neglected. A reproduction of HF results within $\approx 200$ KeV/nucleon at the driplines is obtained, and simple expressions are given for the surface energy and the surface diffuseness parameter. In particular we show that both linear and quadratic terms in $\delta$ and $\Delta R$ are needed to correctly explain the surface term. Moreover, within this analytical mass formula, we show that the neutron skin is essentially determined by the slope of the symmetry energy at saturation, thus confirming earlier numerical results from different groups [36]. Conversely, the surface symmetry energy is shown to be due to a complex interplay of all different local and non-local terms in the energy functional. This implies that constraints on the symmetry energy parameters ($J_{sym}$, $L_{sym}$, $K_{sym}$) from mass measurements might be model dependent and misleading. As a further development of this work, we plan to extend the mass formula to the case of neutron-rich nuclei beyond the dripline in equilibrium with a neutron (and possibly proton) gas. Such a parametrization will allow including modifications of the nuclear surface energy due the presence of continuum
states in nuclear statistical equilibrium models, currently used for different astrophysical applications in supernova and neutron star physics [21].

Acknowledgments

This work has been partially funded by the SN2NS project ANR-10-BLAN-0503 and it has been supported by New-Compstar, COST Action MP1304.

Appendix A. The Gaussian approximation to the isovector energy

In this appendix we develop the 3-dimensional Gaussian $G(r)$ integral

$$E_G = 4\pi \int_0^\infty dr r^2 \mathcal{G}(r) = 4\pi \int_0^\infty dr r^2 A e^{-\frac{(r-R_M)^2}{2\sigma^2}} \quad (A.1)$$

as one-dimensional integrals in order to obtain the isovector surface energy as a function of the nucleus mass and of the effective interaction parameters. As for the symmetric energy, we make the variable change $x = r - R_M$:

$$E_G = 4\pi A \left\{ \int_{-\infty}^{+\infty} dx (x + R_M)^2 e^{-\frac{x^2}{2\sigma^2}} - \int_{-\infty}^{-R_M} dx (x + R_M)^2 e^{-\frac{x^2}{2\sigma^2}} \right\}. \quad (A.2)$$

Since we are interested in the surface energy, we assume that the Gaussian $G(r)$ is zero at the center of the nucleus, such that the second integral in equation (A.2) is negligible with an accuracy $\sim \exp(-R_M^2/(2\sigma^2))$. Then integrating the Gaussian moments straightforwardly lead to

$$E_G = 2(2\pi)^{3/2} \sigma A (R_M^2 + \sigma^2), \quad (A.3)$$

where we have used the expression of the variance as a function of the energy density second derivatives (see section 3). To have $E_G$ as a function of the mass, we just need to express the Gaussian maximum position $R_M$ as a function of A. If we assume $R_M = R$, it reads, using equation (8):

$$E_G = 2(2\pi)^{3/2} \sigma A r_{\text{sat}}^2 \left[ A^{2/3} + \frac{\sigma^2}{r_{\text{sat}}^2} - \frac{2\pi^2}{3} \left(\frac{a}{r_{\text{sat}}}\right)^2 + O\left(\left(\frac{a}{r_0}\right)^4 A^{-2/3}\right) \right] \quad (A.4)$$

In the general case, if we define $\Delta R = R_M - R$, we find additional terms, especially curvature:

$$E_G = 2(2\pi)^{3/2} \sigma A r_{\text{sat}}^2 \left[ A^{2/3} + 2 \frac{\Delta R}{r_{\text{sat}}} A^{1/3} + \frac{\sigma^2}{r_{\text{sat}}^2} - \frac{2\pi^2}{3} \left(\frac{a}{r_{\text{sat}}}\right)^2 + \left(\frac{\Delta R}{r_{\text{sat}}}\right)^2 - \frac{2\pi^2}{3} \frac{\Delta R}{r_{\text{sat}}} \left(\frac{a}{r_0}\right)^2 A^{-1/3} + O\left(\left(\frac{a}{r_0}\right)^4 A^{-2/3}\right) \right] \quad (A.5)$$

Appendix B. The no-skin approximation to the isovector energy

The coefficients entering the isovector energy equations (23) in the no-skin approximation are given by:
$$ C_{\text{surf}}^{\text{IV}} = 3 \left\{ C_{\text{kin}} \left[ \frac{5}{3} \eta_{5/3}^{(0)} \left( \frac{m}{3 m_{\text{sat}}^{*}} + \Delta m_{\text{sat},3} \right) - \left( \frac{\delta m_{\text{sat}}}{3} + \Delta m_{\text{sat},3} \right) \right] - D_{0} \rho_{\text{sat}} + D_{3} \rho_{\text{sat}}^{\alpha+1} \eta_{\alpha+2}^{(0)} \right\}, \quad (\text{B.1}) $$

$$ C_{\text{curv}}^{\text{IV}} = 6 \left\{ C_{\text{kin}} \left[ \frac{5}{3} \left( \eta_{5/3}^{(1)} - \frac{\pi^2}{6} \right) \left( \frac{m}{3 m_{\text{sat}}^{*}} + \Delta m_{\text{sat},3} \right) - \eta_{5/3}^{(0)} \left( \frac{\delta m_{\text{sat}}}{3} + \Delta m_{\text{sat},3} \right) \right] + D_{3} \rho_{\text{sat}}^{\alpha+1} \left( \eta_{\alpha+2}^{(1)} - \frac{\pi^2}{6} \right) \right\}, \quad (\text{B.2}) $$

$$ C_{\text{ind}}^{\text{IV}} = 3 \left\{ C_{\text{kin}} \left[ \frac{5}{3} \left( \eta_{5/3}^{(2)} - \frac{2\pi^2}{3} \eta_{5/3}^{(0)} \right) \left( \frac{m}{3 m_{\text{sat}}^{*}} + \Delta m_{\text{sat},3} \right) - \frac{2}{3} \left( 3\eta_{5/3}^{(1)} - \pi^2 \right) \left( \frac{\delta m_{\text{sat}}}{3} + \Delta m_{\text{sat},3} \right) \right] + \frac{\pi^2}{3} D_{0} \rho_{\text{sat}} + D_{3} \rho_{\text{sat}}^{\alpha+1} \left( \eta_{\alpha+2}^{(2)} - \frac{2\pi^2}{3} \eta_{\alpha+2}^{(0)} \right) \right\}, \quad (\text{B.3}) $$

In these equations, $\Delta m_{\text{sat},3} = (m/m_{\text{sat},n}^* - m/m_{\text{sat},p}^*)/(2I) = (\delta m_{\text{sat},n} - \delta m_{\text{sat},p})/(2I)$. Moreover we have noted $V_{\text{so}} = -mC_{\text{so}}^2/\hbar^2$, where $C_{\text{so}}$ is the usual spin-orbit coefficient of the Skyrme functional [13]. The coefficients entering the isoscalar surface energy equations (17)–(19), as well as the coefficients $\eta_{\gamma}^{(k)}$, can be found in Part I, using $\delta m_{\text{sat}} = (m - m_{\text{sat}}^*)/m_{\text{sat}}^* = (\delta m_{\text{sat},n} + \delta m_{\text{sat},p})/2$, $m/m_{\text{sat}}^* = (m/m_{\text{sat},n}^* + m/m_{\text{sat},p}^*)/2$.

## References

[1] Min L et al 2011 Phys. Rev. C 84 014333

[2] Moller P, Nix J R, Myers W D and Swiatecki W J 1995 ADNDT 59 185

[3] Duflo J and Zuker A P 1995 Phys. Rev. C 52 R23

[4] Aboussir Y, Pearson J M, Dutta A K and Tondeur F 1995 At. Data and Nucl. Data Tab. 61 127

[5] Danielewicz P and Lee J 2009 Nucl. Phys. A 818 36

[6] Lee S J and Mekjian A Z 2011 Phys. Rev. C 82 064319

[7] Nikolov N, Schunck N, Nazarewicz W, Bender M and Pei J 2011 Phys. Rev. C 83 034305

[8] Papakonstantinou P, Margueron J, Gulminelli F and Raduta Ad R 2013 Phys. Rev. C 88 045805

[9] Treiner J and Krivine H 1986 Ann. Phys., NY 170 406

[10] Centelles M, Pi M, Viñas X, Garcias F and Barranco M 1990 Nucl. Phys. A 510 397

[11] Onsi M et al 2008 Phys. Rev. C 77 065805

[12] Potekhin A Y, Fantina A F, Chamel N, Pearson J M and Goriely S 2013 Astron. Astrophys. 560 A48

[13] Aymard F, Gulminelli F and Margueron J 2016 J. Phys. G: Nucl. Part. Phys. 43 045105

[14] Centelles M, Estal M D and Viñas X 1998 Nucl. Phys. A 635 193

[15] Warda M, Viñas X, Roca-Maza X and Centelles M 2009 Phys. Rev. C 80 024316

[16] De J N, Samaddar S K and Agrawal B K 2012 Phys. Lett. B 716 361

[17] Aymard F, Gulminelli F and Margueron J 2014 Phys. Rev. C 89 065807

[18] Ducoin C, Margueron J and Providência C 2010 Eur. Phys. Lett. 91 32001

[19] Ducoin C, Margueron J, Providencia C and Vidaña I 2011 Phys. Rev. C 83 045810

[20] Krivine H and Treiner J 1983 Phys. Lett. B 124 127

[21] CompOSE database (http://compose.obspm.fr/).

[22] Baym G, Bethe H A and Pethick C J 1971 Nucl. Phys. A 175 225

[23] Negele J W and Vautherin D 1973 Nucl. Phys. A 207 298

[24] Raduta A A, Gulminelli F and Aymard F 2014 EPJA 50 24
[25] Aymard F 2015 PhD Thesis University of Caen

[26] Danielewicz P 2003 Nucl. Phys. A **727** 233

[27] Chabanat E, Bonche P, Haensel P, Meyer J and Schaeffer R 1998 Nucl. Phys. A **635** 231

[28] NuDat database (www.nndc.bnl.gov/amdc/)

[29] Reinhard P-G, Bender M, Nazarewic W and Vertse T 2006 Phys. Rev. C **73** 014309

[30] Douchin F, Haensel P and Meyer J 2000 Nucl. Phys. A **665** 419

[31] Myers W D and Swiatecki W J 1980 Nucl. Phys. A **336** 267

[32] Myers W D, Swiatecki W J and Wong C S 1985 Nucl. Phys. A **436** 185

[33] Farine M and Pearson J M 1986 Phys. Lett. B **167** 259

[34] Centelles M, Del Estal M and Viñas X 1998 Nucl. Phys. A **635** 193

[35] Moller P, Myers W D, Sagawa H and Yoshida S 2012 Phys. Rev. Lett. **108** 052501

[36] Ducoin C et al 2011 Phys. Rev. C **83** 045810

[37] Reinhard P-G and Flocard H 1995 Nucl. Phys. A **584** 467

[38] van Giai N and Sagawa H 1981 Phys. Lett. B **106** 379

[39] Cao L G, Lombardo U, Shen C W and Giai N V 2006 Phys. Rev. C **73** 014313