Constraining the surface properties of effective Skyrme interactions

The purpose of this study is threefold: first, to identify a scheme for the determination of the surface energy coefficient a_surf that offers the best compromise between robustness, precision, and numerical efficiency; second, to analyze the correlation between values for a_surf and the characteristic energies of the fission barrier of Pu240; and third, to lay out a procedure how the deformation properties of the Skyrme energy density functional (EDF) can be constrained during the parameter fit. There are several frequently used possibilities to define and calculate the surface energy coefficient a_surf of effective interactions. The most direct access is provided by the model system of semi-infinite nuclear matter, but a_surf can also be extracted from the systematics of binding energies of finite nuclei. Calculations can be carried out either self-consistently (HF), which incorporates quantal shell effects, or in one of the semi-classical Extended Thomas-Fermi (ETF) or Modified Thomas-Fermi (MTF) approximations. The surface energy coefficient of 76 parameterizations of the Skyrme EDF have been calculated. Values obtained with the HF, ETF and MTF methods are not identical, but differ by fairly constant systematic offsets. Despite having some drawbacks compared to the other methods studied here, the MTF approach provides sufficiently precise values for a_surf such that it can be used as a very robust constraint on surface properties during a parameter fit at negligible additional cost. Through the construction of a series of eight parameterizations SLy5s1-SLy5s8 of the standard Skyrme EDF with systematically varied a_surf value it is shown how to arrive at a fit with realistic deformation properties. [Note: The abstract has been abbreviated because of length restrictions imposed by the arXiv. See the paper for the full abstract.]


I. INTRODUCTION
Energy density functional (EDF) methods are versatile tools for the study of nuclear structure and dynamics [1]. Once a parameterization of the EDF has been constructed by selecting terms that incorporate the relevant degrees of freedom and by carefully fixing each that allow for the quantification of confidence intervals of model parameters and the estimate of statistical error bars [5][6][7][8][9][10]. The other is the set-up of more general forms of the Skyrme EDF containing higher-order terms with additional parameters [11][12][13][14][15][16][17]. As both developments substantially increase the numerical cost of the parameters' adjustment, it is of advantage to have efficient methods to calculate the data and pseudo-data used during the fit.
In the past, many authors have discussed the crucial role that an accurate adjustment of the surface energy coefficient of an EDF's parameterization plays for the deformation properties of nuclei [7,9,[18][19][20][21][22]. Generally speaking, when increasing a surf , deformation energy surfaces become stiffer, fission barriers higher and the excitation energy of fission isomers and other superdeformed states more elevated. This, in turn can then be used to fine-tune the EDF. First pioneering calculations of fission properties indicated that the two early Skyrme parameterizations SIII [23] and SkM [24] give a fission barrier height for 240 Pu that is too high or too low by roughly a factor of two, respectively. This disagreement clearly exhibited within semi-classical calculations [25] has triggered the adjustment of the SkM * parameterization, for which the parameters of the momentum-dependent part of SkM have been modified and fine-tuned in order to change the height of the semi-classical fission barrier without changing other properties of infinite nuclear matter [18]. Even after more than 30 years, SkM * is a still often used reference for the description of fission phenomena [26][27][28][29][30][31].
The aim of the present study is threefold. First, we will set up and benchmark an efficient and robust method to calculate the surface energy coefficient a surf as defined in Ref. [32] for modern Skyrme EDFs. Second, we will analyze the correlation between the value for a surf and characteristic energies of the deformation energy landscape of 240 Pu. Third, we will outline how a constraint on the value of a surf can be incorporated into the parameter adjustment.
This article is organized as follows. In Section II we present the three methods based on a one-dimensional model for semi-infinite nuclear matter that will be used to evaluate the surface energy coefficient of the Skyrme EDF: the self-consistent Hartree-Fock method (HF), the extended Thomas-Fermi method (ETF), and the Modified Thomas-Fermi method (MTF). Section III recalls the "leptodermous protocol" of Reinhard et al. [33] to extract the surface energy from an analysis of very large fictitious nuclei. We will use this approach with the ETF method. A systematic comparison of the results provided by these four possibilities is presented in Section IV. In Sect. V, we discuss to what extent the such extracted surface energy coefficient is linked to characteristic energies of the fission barrier of 240 Pu. Finally, Section VI presents the adjustment of a series of standard Skyrme parameterizations with systematically varied a surf , called SLy5s1-SLy5s8. A summary of our discussion and per-spectives will be given in Section VII.

II. SURFACE ENERGY IN SEMI-INFINITE NUCLEAR MATTER
A. The Skyrme EDF The starting point for the derivation of a Skyrme EDF is usually an effective two-body interaction with parameters that may depend on the density of the system. It is important to recall that this effective interaction is only used as a generator for the general form of the EDF and allows for writing the coupling constants that weight these scalars as functions of a limited number of parameters. It is then common to disregard certain terms in the functional or to relax some of the interdependences between the coupling constants, which breaks the one-toone-correspondence between the EDF and the underlying effective interaction used to construct it [1,11,17,34].
The present standard form of the Skyrme EDF is motivated by the use of a density-dependent two-body interaction [1], leading to a bilinear EDF with densitydependent coupling constants. Most of the parameterizations discussed below will be of that type. However, we will also compare with results obtained for EDFs derived from interactions with two-, three-and sometimes even four-body terms but density-independent coupling constants [16,17,35]. Such functionals will be referred to as trilinear or quartic EDFs, respectively.
In general, the total energy can be written as the sum of five terms [1]: the kinetic energy E kin , a potential energy functional E Sky that models the strong interaction in the particle-hole channel, a pairing energy functional E pairing , a Coulomb energy functional E Coulomb and a correction term E corr that approximately removes the excitation energy due to spurious motions caused by broken symmetries E = E kin + E Sky + E pairing + E Coulomb + E corr . (1) The model systems that we will use to extract the surface properties are semi-infinite nuclear matter and giant spherical nuclei, for which only the first two parts of the energy are taken into account. As a consequence of the Skyrme interaction being a contact force, the corresponding EDF then can be written in the form of an integral over a local energy density where the Skyrme part E Sky (r) can be further decomposed into central (E C t ), spin-orbit (E LS t ) and tensor (E T t ) terms E Sky (r) = t=0,1 that are either composed entirely of isoscalar densities (t = 0) or that contain bilinear combinations of isovector densities (t = 1). The Skyrme energy functional as such is constructed to be an isoscalar. The physics contained in the Skyrme functional has been discussed in great detail in the literature [1,11,17,34,36], and here we will use standard notations for local densities and coefficients of the functional. Its complete form discussed in these papers contains many terms, only a small subset of which are present for the systems we consider here. In particular, all time-odd densities are zero because of the time-reversal invariance we impose, and many components of the time-even vector and tensor densities are zero because of spatial symmetries.

B. Model of semi-infinite nuclear matter
The surface energy is often extracted from an idealized one-dimensional model of semi-infinite nuclear matter originally developed by Swiatecki [37] and revisited in [32]. One considers a medium where the local densities are constant along the x and y directions, but vary along the z direction. The corresponding profile of the local proton and neutron matter densities will be noted as ρ q (z), q = n, p. Deep inside the matter for z → −∞, one expects that ρ 0 (z) → ρ sat , i.e. the equilibrium density of infinite nuclear matter, and that E(z) → a v , i.e. the volume energy per particle at saturation. Far outside the matter, i.e. for z → +∞, one has ρ 0 (z) → 0. For sake of compact notation, we will drop the z-dependence of the densities from hereon whenever possible.
Within this one-dimensional model, the only nonvanishing components of the cartesian spin-current tensor density J 0,µν are J 0,xy = −J 0,yx . As a consequence, only the J 0,z component of the vector part J 0 of the full spin-current tensor as defined in [11,36] is non-zero.
We will here focus on the discussion of the surface energy coefficient a surf that is related to the surface energy of symmetric semi-infinite nuclear matter. In this system, proton and neuton densities are equal, i.e. ρ n (z) = ρ p (z) = 1 2 ρ 0 (z) and similar for the other densities. As a consequence, only those terms in the EDF that are entirely composed of isoscalar densities have to be considered. This enormously simplifies the energy density (2) that can be reduced to 1 1 Note that no unique definition of the coupling constant C J t of the tensor terms can be found in the literature. We use here the convention of Ref. [11]. Others might differ by a factor two.
where the density-dependent isoscalar effective mass is given by the ratio The expression for the Skyrme EDF provided by Eq. (5) covers a vast number of different parameterizations that have been lately used in the literature. We use a notation where the coupling constants for the bilinear, trilinear and quartic terms are denoted with C 0 , B 0 and D 0 , respectively.
In the most-widely used standard form of the Skyrme EDF only the bilinear terms are considered, i.e. all B 0 = D 0 = 0, and the coupling constant C ρ 0 is made explicitly density dependent by multiplying it with [1 + c ρ α 0 (r)], where the parameter c controls the relative weight of the density-dependent part of the coupling constant.
The J 2 0 terms bilinear in the spin-current density will be called "tensor terms" in what follows. For a majority of the widely-used parameterizations the coupling constant C J 0 of the tensor terms is set to zero, such that E T t = 0 in Eq. (4). Many past semi-classical and HF calculations of semi-infinite nuclear matter, however, have neglected these terms also for those parameterizations for which they are to be taken into account.
The possibility of replacing the density dependence of C ρ 0 by trilinear (B 0 = 0) and quartic (D 0 = 0) terms, and where all coupling constants are derived from an underlying Skyrme 2+3+4-body Hamiltonian, has been considered recently with the goal of constructing well-defined EDFs for use in beyond-mean-field methods [16,17]. This extended form of the Skyrme EDF will also be considered in the present work.
The assumptions made when setting-up the Skyrme EDF have been discussed in great detail in the literature [1,4,11,17,34,36] and will not be recapitulated here. Neither will we repeat the discussion concerning the different possible definitions for the surface energy. For a detailed discussion, we refer to the original articles of Myers and Swiatecki [38,39] and the more recent ones by Pearson et al. [40,41], Brack et al. [25], Kolehmainen et al. [42], Centelles et al. [43] or Douchin et al. [44].

C. Hartree-Fock calculations
The first one of the methods we use to calculate semiinfinite nuclear matter is the self-consistent mean-field approximation, usually called Hartree-Fock (HF), with a treatment along the lines of Refs. [32,40]. In this case, one considers the quantity [32] which represents the energy per unit of surface for a piece of semi-infinite nuclear matter described by a density which is constant in the x and y directions and extends from −L to +L in the z direction with the conditions ρ 0 (−L) = ρ sat and ρ 0 (L) = 0 for L → +∞. After solving the mean-field equations, the surface energy coefficient, denoted as a HF surf , can be extracted using with the parameter r 0 being defined through the condition 4 3 πr 3 0 ρ sat = 1. As a particular feature of such quantal calculation, one observes so-called "Friedel oscillations" [45,46] of the density ρ 0 (z) in the vicinity of the surface inside the matter. Since these oscillations are only very slowly damped, a reliable calculation of a HF surf requires quite large an interval in z direction.
Our HF values for a surf often differ slightly from those given by other groups in the past [18,[47][48][49][50][51], see the supplementary material [52]. On the one hand, this underlines the numerical difficulties of determining a precise value for a surf . On the other hand, as said before, the contribution from the J 2 tensor terms, which are present for a subset of the parameterizations, has been omitted in most of the earlier published work.

D. Semi-classical calculations
As a second method to determine a surf we use the semiclassical Extended Thomas-Fermi (ETF) approach up to order 4 [25]. Values obtained with this method will be denoted as a ETF surf in what follows. In this semi-classical framework, the local densities ρ q (z) are modeled by a three-parameter modified Fermi function. The kinetic and spin-current densities entering the total energy density E(z) are obtained from an expansion of the so-called single-particle Bloch density matrix in powers of around its Thomas-Fermi value as originally proposed by Wigner and Kirkwood. In short, the Bloch density matrix is the coordinate-space representation of the statistical operator e βĥ constructed from the HF single-particle Hamiltonianĥ expressed in its eigenbasis [53], and related to the usual coordinate space representation of the density matrix ρ(r, r ′ ) through a Laplace transform [25,53]. Ultimately, this leads to expressions for τ 0 and J 0 as functions of the local density ρ 0 and its derivatives, where it is customary to separate the contributions of different (even) power in in the expansion J 0 = J [2] 0 + J [4] 0 .
The complete expressions for τ [0] 0 (which is simply the kinetic energy of a non-interacting Fermi gas), τ [2] 0 , τ [4] 0 , J [2] 0 and J [4] 0 have been given by Brack et al. [25] and more recently by Bartel et al. [54]. Note that, in the case where tensor terms are included in the Skyrme EDF, the standard expressions given in the early articles have to be modified taking into account the results of Bartel et al. [55] for the contribution of the J 2 terms. Minimizing the surface energy as calculated with an expression equal to the one given by Eq. (8), one obtains the parameters of the assumed Fermi-type density profiles from which the surface energy coefficients can then be calculated.

E. Modified Thomas-Fermi approximation
The ETF approximation provides expressions for the kinetic and spin-current densities in terms of the nucleon density and its derivatives. The Modified Thomas-Fermi (MTF) approximation developed by Krivine and Treiner [56] consists in using an ETF expansion limited to order 2 where the coefficients are modified to simulate the order 4 as well as the 2 effective mass contributions. The MTF form of the kinetic and spin-current densities can then be written as τ [2] with the coefficients α = 3 5 , β = 1 18 and γ = 1 3 , chosen for a resonable reproduction of the total energy of finite nuclei [56]. In local density approximation, the Fermi momentum is given by k F = ( 3 2 π 2 ρ 0 ) 1/3 and the spinorbit field W 0 is defined as usual as W 0 = ∂E/∂J 0 .
The coefficients α and γ of the MTF expressions thus keep their values obtained from ETF, whereas the value of β is changed from 1/36 (obtained within ETF) to 1/18 (used in the MTF method).
The interest of the MTF method is that in the case of symmetric semi-infinite matter, the approximations provided by Eqs. (11), (12) and (13) lead to an analytically solvable Euler-Lagrange equation for the density profile without the need for carrying out a variational calculation numerically. It is noteworthy that the entire density profile is varied in the MTF method, whereas in the ETF approach only the parameters of a predefined Fermi function are optimized to give the lowest binding energy.
Following the derivation outlined in Ref. [57], the semiclassical energy density of symmetric matter can then be simply written as the sum of two terms where Compared to what is found elsewhere in the literature, the expressions (15) and (16) have been complemented in order to take into account the possible 3-and 4-body contributions to the Skyrme EDF (5) as introduced in Refs. [16,17]. In addition, in order to take into account the tensor terms in the Skyrme EDF, the vector spincurrent density is obtained through the protocol of Bartel et al. [55], which leads to a redefinition of the V so [ρ] term that, compared to the original work by Treiner and Krivine [57], contains additional terms.
Inserting the analytical solution for the density profile ρ 0 (z) into the expression for the surface energy, one obtains after some further analytical manipulations a compact expression for the surface energy coefficient (for details see [57]), denoted a MTF surf from hereon, As the MTF method differs from ETF approach only by the values of the coefficients, this expression has exactly the same form as the Wilets formula [58,59] derived much earlier.
Hereafter, we will call Eq. (22) the MTF pocket formula in order to underline the simplicity of its use to calculate a surf compared to the HF and ETF methods described above. In its original MTF form without the contribution from the tensor terms, this expression has been occasionally used to analyze the surface properties of Skyrme parameterizations [57,60,61].

III. SURFACE ENERGY FROM SEMICLASSICAL CALCULATIONS OF LARGE NUCLEI
In Ref. [33], Reinhard et al. have proposed a protocol to extract the nuclear liquid-drop coefficients that correspond to an EDF parameterization from a leptodermous expansion based on mean-field calculations of very large fictitious spherical nuclei calculated without Coulomb interaction and pairing correlations. For the determination of the (isoscalar) surface energy coefficient, the calculations can be limited to symmetric nuclei with N = Z = A/2. In this case, the calculated binding energy of the finite nuclei E(A) can then be developed into i.e. a volume, surface, and curvature term plus another one that is proportional to A 0 = 1. The latter has not been considered in Ref. [33]. Despite its unrealistic limit for A → 0, such term emerges naturally in the leptodermous expansion of the nuclear binding energy as a second-order correction to the curvature energy [62,63], or when replacing the geometric surface energy of the standard liquid-drop model by a double-folding integral [64,65] as it is done in the finite-range liquid-drop and droplet models (FRLDM) [66].
For consistency with semi-infinite matter calculations, the finite nuclei are calculated without center-of-mass correction, irrespective of the scheme used during the fit of the parameterization used.
In order to extract the surface energy coefficient a surf of Eq. (23), one first defines an effective surface energy coefficient a eff surf (A) of a given nucleus of mass A by reshuffling the expansion (23) The volume energy coefficient a vol is provided by the energy per particle E/A of symmetric infinite nuclear matter at saturation density. By fitting a second-order polynomial in A −1/3 to the calculated values for a eff surf (A), one obtains the coefficients a surf , a curv and a 0 .
In order to disentangle the surface energy unambiguously from the higher-order terms in the liquid-drop formula (23), we calculate 20 nuclei with very large mass numbers in between 1 200 ≤ A ≤ 200 000, similar to what was done in Ref. [33]. There, however, nuclei were calculated self-consistently, which required to remove shell effects by subtracting the shell correction as obtained from the self-consistent single-particle spectrum from the total binding energy. Here, we use the semi-classical ETF approach to calculate binding energies instead, such that there is no need to eliminate shell effects. The expansion in terms of powers of A −1/3 (or inverse nuclear radius) is thus more stable and the LDM parameters thus more precisely extracted when extrapolating to (semi-)infinite nuclear matter in the limit A −1/3 → 0. The such determined values for a surf , however, should be compared to the ETF results obtained for semi-infinite-matter, and not the HF results as the ones extracted in Ref. [33].

IV. DETERMINATION OF a surf
In this section we present a systematic comparison of results for the surface energy coefficient obtained with the four aforementioned methods. One of our main goals is to check whether the computationally friendly MTF pocket formula provides a reliable estimate for the surface energy coefficient as extracted from the theoretically more advanced, but numerically more involved, HF or ETF methods. For that purpose, we use a large set of Skyrme parameterizations.

A. Parameterizations considered
There is a large number of Skyrme parameterizations that can be found in the literature. Only few of them, however, are frequently used in production runs. The simplicity and popularity of the Skyrme EDF has led to quite large a number of "experimental" fits of parameterizations that were carried out for one or the other very specific study. In particular, there are many "families" of fits which explore the influence of variations of details of the parameterizations on their predictive power. We will profit here from this large number of parameterizations as it allows us to cover large intervals of values for infinite nuclear matter properties, which in turn can reveal possible correlations of differences between the four methods to determine a surf and other global features of the parameterizations.
Here we give a brief overview over the main features and particularities of the parameterizations considered here (which are all of the standard density-dependent bilinear form unless specified otherwise), • SIII [23]: the coupling constants C ρ t of this early, but still sometimes used, parameterization are linear functions of the density ρ 0 ; • Ska [67], SGI, SGII [68], SkM [24], SkM* [18]: these are examples of early standard parameterizations with a density dependence with α < 1, taking values of either 1/3 or 1/6 as almost all parameterizations listed below, unless otherwise specified; • The SLy family of fits: more recent and widelyused examples of standard parameterizations are the Saclay-Lyon fits SLy4-SLy7 [69] that differ in options for center-of-mass correction and tensor terms, SLy5* [70], a recent refit that suppresses the finite-size spin-instability of the original SLy5 parametrization, and SLy4d [71], a refit of SLy4 without any center-of-mass correction built for the purpose of TDHF calculations; • The TIJ family of fits [11]: for these parameterizations, also fitted within the Saclay-Lyon protocol, the coupling constants of the tensor terms were systematically varied over large intervals; • SLy5+T [72], SLy4T, SLy4T min [73], SLy4T self , TZA [74]: these are further fits based on the Saclay-Lyon family of fits with either added or modified tensor terms; • f 0 and f ± [75]: these constitute another series of variants of Saclay-Lyon parameterizations that explore different values for the splitting of the effective mass of protons and neutrons in asymmetric matter. Unlike the other standard parameterizations considered here, their coupling constants C ρ t have two density dependences with powers α 1 = 1/3 and α 2 = 2/3, respectively. This allows for the decoupling of nuclear matter properties that cannot be chosen independently for standard parameterizations with just one density dependence [76]; • The SLyIII.x family of fits [77]: these are yet another series of variants of Saclay-Lyon parameterizations that were built specifically for the purpose of regularized beyond-mean-field calculations. All have the same linear density dependence of the C ρ t coupling constant as SIII. Their isoscalar effective mass m * 0 /m = x has been constrained to the values 0.7, 0.8, 0.9 and 1.0 in units of the nucleon mass; • SaMi [78]: this recent parameterization has been adjusted within a SLy-inspired protocol with the aim of improved spin-isospin properties; • SkI3, SkI4 [79], SkO, SkO' [80]: these fits are examples of standard parameterizations with a generalized isospin dependence of the spin-orbit part of the EDF; • The SV family of fits [81]: this series of parameterizations was constructed for the purpose of studying the correlation of nuclear matter properties with other observables. Starting from the reference parameterization called SV-bas, nine others were constructed by varying the incompressibility K ∞ , isoscalar effective mass m * 0 /m, symmetry energy J and the sum-rule enhancement factor κ v one by one while keeping the others constant. A best fit, called SV-min, is also considered; • The UNEDF family of fits: these parameterizations differ by the selection of data considered in the fit protocol. Compared to UNEDF0 [6], for example, UNEDF1 [7] and UNEDF2 [9] are also adjusted to reproduce the excitation energy of the fission isomer of 240 Pu. We also consider the parameterization UNEDF1 SO [82] that corresponds to a readjustment of the spin-orbit coupling constants of UNEDF1, thereby improving the spectroscopy of very heavy nuclei at the expense of a much lower value for a surf ; • The BSk family of fits: we also consider several representative parameterizations from the series of large-scale Bruxelles-Montreal Skyrme-HFB mass fits that stay within the standard form of the Skyrme EDF with only the coupling constants C ρ t being density dependent [21,49,[83][84][85][86][87][88], among which BSk14 [21] has also been fitted to fission barriers; • LNS [89], NRAPR [90], NRAPRii [91]: in one way or the other, these parameterizations of the standard Skyrme EDF were adjusted to reproduce nuclear matter properties as predicted by ab-initio methods. While LNS was adjusted to reproduce a large variety of nuclear matter results from a Brueckner-HF calculation, NRAPR has been fitted to the density dependence of the energy per particle as obtained from ab-initio calculations of nuclear and neutron matter. We also consider the parameterization NRAPRii with doubled strength W 0 of the spin-orbit interaction compared to NRAPR, as suggested by the authors of Ref. [91]; • KDE0v1 [61]: this standard parameterization has been adjusted to reproduce a large number of empirical nuclear matter data; • SQMC700 [92]: this parameterization of the standard Skyrme EDF has been derived as the nonrelativistic mean-field limit of a quark-mesoncoupling model [92]; • S1Sd, S1Se [93]: the parameters of these two standard density-dependent Skyrme parameterizations have been adjusted to reproduce total binding energies of doubly-magic nuclei as predicted by the Gogny force D1S. For S1Sd the tensor terms were included, whereas for S1Se they were neglected; • S3Ly family of fits [35]: we also include a few representative examples from the series of fits of extended Skyrme EDFs that add central threebody terms with gradients to a standard densitydependent two-body Skyrme EDF, and which were carried out within a modified Saclay-Lyon fit protocol. These fits do systematically cover a wide range of values for the isoscalar effective mass m * 0 /m and incompressibility K ∞ . For example, S3Ly71260 is a parameterization with m * /m 0 = 0.71 and K ∞ = 260 MeV; • SLyMR0 [16], SLyMR1 [94]: we also include two recent parameterizations built for the purpose of spuriousity-free beyond-mean-field calculations: SLyMR0, which combines the non-densitydependent part of the standard 2-body central and spin-orbit Skyrme interaction with gradient-less 3and 4-body terms, and SLyMR1, where the 4-body terms are replaced by the 3-body terms with gradients as introduced in Ref. [17].
The main interest of the four rarely-used parameterizations KDE0v1, LNS, NRAPRii and SQMC700 is that they were recently shown to be consistent with a large set of pseudo-data for infinite symmetric nuclear matter [95]. Their predictive power for finite nuclei, however, is rather limited [91]. Most, if not all, other parameterizations of the standard density-dependent Skyrme EDF listed above provide a much better description of finite nuclei than these, but in turn are incompatible with some of the presently accepted values for the empirical properties of nuclear matter [95]. It should be stressed, however, that none of the nuclear matter properties is a real observable, as in one way or the other they all have to be extracted in a model-dependent way from data.
B. Surface energy coefficients Figure 1 shows the surface energy coefficients a surf obtained from the HF, ETF and MTF methods for the list of parameterizations given above. Two sets of ETF values, one obtained from the calculation of semi-infinite nuclear matter (open blue squares), the other extracted from calculations of large finite nuclei (blue dots), are shown.
The first observation that can be made is that the values of the surface energy coefficient spread over a relatively large interval, from about 15.5 MeV to about 19.5 MeV for the results given by the HF calculations. As many parameterizations were obtained within dissimilar protocols, it is difficult to correlate the value for a surf with other properties of the respective parameterizations. A few correlations that can be unambiguously identified are that a surf depends on the presence and size of tensor terms (compare SLy4, SLy5, SLy5T, and the TIJ), the size of the spin-orbit term (compare NRAPR and NRAPRii), the scheme for center-of-mass correction used (as already pointed out in [20]; compare SLy4, SLy6 and SLy4d). One can expect that there are further correlations between a parameterization's value for a surf and its other properties. The analysis of their origin and nature, i.e. if they are a physical necessity or rather an accidental consequence of either a specific fit protocol or an overconstrained form of the EDF, however, is beyond the scope of the present study.
A second observation is that the differences between the values for a HF surf and a ETF surf on the one hand, and the differences between a HF surf and a MTF surf on the other hand, are fairly constant and almost independent on the nature of the EDF considered (density-dependent bilinear, trilinear, or even quartic) and the properties it provides for infinite nuclear matter.
A third observation is that for all parameterizations the discrepancy between the two different ETF values is so small that it can be hardly resolved on the Figure, which confirms that the leptodermous protocol of Ref. [33] offers a reliable alternative to calculations of semi-infinite nuclear matter, with the remaining differences being on the order of 10 keV.
The size and parametrization-dependence of the little remaining scatter between the three methods to determine a surf from semi-infinite nuclear matter calculations can be better resolved on Fig. 2 TZA   T66   T64   T62   T46   T44   T42   T26   T24   T22   17   In any event, Fig. 2 indicates that, for the purpose of calculating a surf , MTF is almost as good an approximation to HF as ETF. However, the two semi-classical methods differ among each other on a much larger scale as one systematically overestimates the HF value, whereas the other systematically underestimates it.
A closer examination of the parameterizations for which the scatter is largest indicates that ∆a HF-MTF surf may strongly depend on the isoscalar effective mass, see for example the SLyIII.x and BSk series for which m * 0 /m ranges between 0.7 to 1.05).
In some cases one also finds small differences between parameterizations with similar effective mass that are correlated to the strength of the tensor terms. This is exemplified by the comparison of SLy4, SLy5, SLy5+T and the TIJ series that were adjusted within almost the same protocol, have practically the same effective mass, but differ in the size and sign of the tensor coupling constants C J t . For these cases the difference between ETF and MTF seems not to be much affected, but that the differences between HF and ETF and also between HF and MTF are not the same. In general, Fig. 2 suggests that without tensor terms, ETF and MTF have the same offset from HF, but in opposite directions. In the presence of tensor terms the ETF values for a surf come closer to the ones from HF, whereas MTF moves further away. This finding points to limitations of modeling the spincurrent density in semi-classical methods.
Comparing SLy5 and SLyIII.0.7, which have similar effective mass and tensor terms, indicates that ∆a HF-MTF surf might in addition also depend on the power of the density dependence of C ρ 0 [ρ 0 ]. For α = 1 (SLyIII.0.7), MTF values for a surf are much closer to the ones from HF than for α = 1/6 (SLy5). We will come back to this below.
One can expect that there might be further weak correlations between a parameterization's value for ∆a HF-MTF surf and its other properties, but these cannot easily be identified even within this large set of parameterizations.
In order to analyze further the correlation between ∆a HF-MTF surf and nuclear matter properties, parameterizations with systematically varied nuclear matter properties are needed. Such fits are provided by the series of SV parameterizations by Klüpfel et al. [81]. During their adjustment, the incompressibility K ∞ , isoscalar effective mass m * 0 /m, symmetry energy coefficient J and the Thomas-Reiche-Kuhn enhancement factor κ v have been separately varied while keeping the other properties constant.
The left panel of Fig. 3 shows the corresponding values for a surf obtained with the four schemes to calculate the surface energy coefficient introduced above. Two correlations for the size of a sym itself, independent of the method of its determination, become immediately obvious: a surf increases with increasing incompressibility K ∞ and also with increasing symmetry energy coefficient J. Both, however, might be particular to the specific fit protocol of this series of parameterizations.
More importantly, for the series with varied effective mass m * 0 /m, the different methods to determine a surf give a different trend, as already indicated by the analysis of the large set of parameterizations in Fig. 2. While the HF and ETF values are fairly independent on m * 0 /m, the MTF values rapidly increase with decreasing effective mass. The method dependence of the a surf values becomes more obvious in Fig. 3, where are plotted the differences between the methods in the same way as in Fig. 2.
The right panel of Fig. 3 also reveals that there is a slight dependence of both ∆a HF-ETF surf and ∆a HF-MTF surf on the incompressibility K ∞ , which to a large extent explains the difference of ∆a HF-MTF surf between SLy5 and SLyIII.0.7 found in Fig. 2 and attributed to the power of the densitydependent term above. The values for K ∞ of these two parameterizations are quite different, K ∞ = 230 MeV for SLy5 and K ∞ = 361.3 MeV for SLyIII.0.7, which is a consequence of the correlation between m * 0 /m, K ∞ and the power of the density-dependence of the C ρ 0 coupling constant analyzed for example in Ref. [76]. The K ∞ value of SLyIII.0.7 is far outside the range covered by the SV parameterizations of Fig. 3 (and therefore also far from the empirical value, which is unavoidable for standard parameterizations with linear density dependence [76]). Going from SV-K218 to SV-K241, one finds that ∆a HF-MTF surf /∆K ∞ ≈ 0.002. Assuming the same slope when going from SLy5 to SLyIII.0.7, then their ∆a HF-MTF surf should differ by about 260 keV, which is indeed the case.
Altogether, we find a reasonably parameterizationindependent behavior of the three models with a quantitative difference that is almost constant as long as K ∞ and m * 0 /m are constrained to a small interval, which is usually the case within a given fit protocol. From this, one can conclude that the Wilets pocket formula from Eq. (22) can be safely used in a fit protocol to determine the coupling constants of future Skyrme EDFs.

C. Density profiles in semi-infinite nuclear matter
To examine the origin of the systematic differences found above between the three models used to calculate semi-infinite nuclear matter, Fig. 4 displays the radial profiles of the mass density ρ(z), kinetic density τ (z), and the z-component of the spin-current density J z (z) as obtained from each method. For each model, the profiles are positioned such that the lower boundary −L ′ of a piece of the surface with some sufficiently large, but otherwise arbitrarily chosen, particle number A calculated as coincides for all of the three models. 2 The parameter r 0 is the same as in Eq. (8). With this, the density profiles can be compared exactly as those of finite nuclei with same particle number. The origin z = 0, however, has been chosen to correspond to the position of the sharp surface of a piece of saturated nuclear matter with constant density inside that is placed in precisely the same manner as the three density profiles shown in Fig. 4, i.e. −L ′ = A/(4πr 2 0 ρ sat ) = A/[3 1/3 (4πρ sat ) 2/3 ]. As for finite nuclei, the density and kinetic density are bulk properties that approach a saturation value inside the matter and fall off at the surface. By contrast, the spin-current density is peaked on the surface and approaches zero inside and outside the matter. Indeed, J z is exactly zero in infinite homogeneous nuclear matter when calculated in mean-field approximation.
The amplitude of the Friedel oscillations exhibited by all three HF densities is comparatively small and barely visible when plotting the entire density profile. However, as indicated by the inserts, the oscillations with a wavelength of about 2.4 fm are only slowly damped and reach far inside.
Each of the three approaches to calculate semi-infinite matter provides a slightly different profile of the mass density distribution ρ(z) at the surface. They differ in the steepness, diffuseness, and also in the position of the inflection point. This is not very surprising, as in one way or the other the ETF and MTF methods are constructed to provide an approximation to the total energy of an HF calculation, and not to reproduce its densities. While in the ETF method only the parameters of a Fermi function parameterizing the profile of the local matter density ρ(z) are variationally optimized, the entire density profile is implicitly varied in the MTF approach.
The ETF density ρ(z) follows quite closely the mean of the oscillations of the HF density, such that the two are difficult to distinguish in Fig. 4. By contrast, the profile of the MTF density is visibly different from the other two by being much more asymmetric around the surface. Inside the matter, the MTF density approaches the saturation value slower than the ETF density or the mean of the oscillating HF density, whereas outside the matter it falls much quicker to zero than the other two. The same behavior has already been found for finite nuclei in the seminal papers on the MTF method [56,57]. This asymmetry has some consequences for the corresponding kinetic and spin-current densities. On the one hand, it generates a visible bump in the tail of τ (z), and also leads to a very asymmetric shape of J z (z).
A similar, but much less pronounced, bump is also found in the tail of the ETF kinetic density τ (z). In this case, it is generated by the terms of order 4 in the semi-classical expansion, which are also at the origin of the bump that the ETF spin-current density J z (z) exhibits at the same position. The appearance and size of the latter is parameterization-dependent, and can be correlated to the relative sign and size of spin-orbit and tensor terms.
While these differences in the density profiles easily explain why the three methods deliver slightly different results for the surface energy coefficient, it is more difficult to correlate the dissimilarities of the density profiles with the systematic differences found for the a surf values. In any event, it has to be recalled that the ETF and MTF approaches are approximations to the HF calculation of the total energy of the nuclear system. The surface energy coefficient defined through Eq. (8), however, is the difference between two energies that are typically two orders of magnitude larger, such that the differences found between the various schemes to calculate them is beyond the third significant digit, which should not be unexpected.

V. RELATION BETWEEN a surf AND DEFORMATION ENERGY SURFACES
It is well understood that the characteristic features of the deformation energy surfaces of heavy nuclei, i.e. the excitation energies of secondary minima and the height of barriers separating the minima and stabilizing the nucleus against fission, are strongly correlated with the surface properties of the effective interaction [18,22]. An emblematic example is the "double-humped" fission barrier of 240 Pu [96] that we will use here as illustrative example. As a rule of thumb, the larger the surface energy coefficient a surf , the higher the excitation energies and barriers.
In the liquid-drop model, the surface energy of a nucleus is simply provided by the product of a surf , the size of the nucleus' surface and a universal factor. Assuming that the deformation of the fission isomer and the top of the barriers turn out to be at the same deformation for all parameterizations (such that the nuclear surface is of comparable size), one would then naively expect that their energy is linearly correlated with the surface energy coefficient. In addition, the larger the deformation, the steeper should be the slope.
It has to be recalled, however, that the characteristic features of the energy surfaces, in particular the groundstate deformation and the presence of secondary minima, are caused by shell effects [96,97]. As a consequence, the variation of shell effects with deformation is as important for the observable excitation energies as is the smooth variation of the of the liquid drop surface energy with deformation, which on its own would just give all actinide nuclei a spherical ground state and one structureless broad fission barrier. And indeed, the fact that the amplitude of the variation of shell effects with deformation that leads to the characteristic double-humped fission barrier of most actinides is correlated with the effective mass and spin-orbit strength has been pointed out already very early by Tondeur [98].
For a subset of the parameterizations employed above, we have carried out calculations of the complete static fission barrier of 240 Pu, which is an often used benchmark for such studies [1, 18-20, 28, 99-102]. In order to avoid a readjustment of the pairing strength and the ambiguities related to it, we have limited the analysis to parameterizations of similar effective mass. Most were adjusted with a variant of the Saclay-Lyon protocol. In the figures we will only distinguish between series of fits, which are the SLyx family, (SLy4-7, SLy5*, SLy4d), the fx family (f 0 , f ± ), the TIJ family (T22, T24, T26, T42,  T44, T46, T62, T64), and the SLyxT family (SLy5+T, SLy4T, SLy4T min , SLy4T self , TZA). In addition, we have used the classic SkM and SkM* parameterizations. As it turns out, this subset is sufficient for a conclusive analysis.
In all cases, we use "surface pairing" with strength −1250 MeV fm −3 for protons and neutrons, and a soft cutoff at ±5 MeV above and below the Fermi energy as defined in Ref. [103].
Calculations are carried out with the most recent versions of the Ev8 [34] and Ev4 [104] codes. Both use the same 3D coordinate-space representation, where the single-particle wave functions are discretized on an equidistant mesh in a box. The two codes differ by the symmetries they impose on the nuclear shapes. Ev8 assumes three plane reflection symmetries, which reduces the calculation to 1/8 of the full box, but is still sufficient to describe triaxial shapes. By contrast, Ev4 assumes only two plane reflection symmetries, which then also permits to describe (not necessarily axial) octupole deformed shapes. The Ev8 box has n x × n y × n z = 20 × 20 × 30 points of distance 0.8 fm. The Poisson equation for the Coulomb field is solved in a 50 × 50 × 50 box for improved precision of the Coulomb energy at large elongation. The Ev4 box doubles the z dimension. With these choices, the deformation energy reaches a precision of about 100 keV independent on the nuclear shape [105].
Up to the fission isomer, the shapes along the static fission path are reflection symmetric, whereas beyond the fission isomer the shapes quickly become massasymmetric with increasing quadrupole deformation. The inner barrier is triaxial, everywhere else shapes along the path are axial. Figure 5 shows the excitation energy of the superdeformed fission isomer and the height of the inner and outer barrier as a function of the surface energy coefficient calculated in three different ways.
For the experimental energies, which are indicated by grey horizontal bars in Fig. 5 [108], the first barrier has 5.8 MeV and the second barrier 5.3 MeV. A more recent paper from the same group [49] quotes < 5.7 MeV for the inner and 4.5 MeV for the outer barrier. In their paper on the fit of UNEDF1, Kortelainen et al. [7] quote the values 6.05 MeV for the inner barrier and 5.15 MeV for the outer one, as recommended by the RIPL-3 database [109,110]. The scatter of the experimental data, however, is much smaller than the scatter in the theoretical results of Fig. 5. In any event, most parameterizations largely overestimate the experimental values, which is the principal motivation for the current efforts to improve the EDFs in that respect.
Note that none of the parameterizations does simultaneously describe all of the three properties. The three parameterizations that give a reasonable description of E iso of about 2.7 MeV (and which are SkM*, SLy6 and SLy7), still overestimate the inner and outer barrier heights by more than 1 MeV.
All deformation energies increase with the value of the surface energy coefficient a surf as expected. However, the correlations are not strictly linear. Instead, there is a large scatter around the global trends, in particular for the height of the inner barrier. This is an immediate consequence of the shell effects not being the same for all parameterizations. Indeed, it has been demonstrated in Ref. [22] that the contribution of the shell correction as deduced from self-consistent mean-field calculations of the excitation energy of the fission isomer of nuclei in the actinide region, including 240 Pu, can vary by more than one MeV when going from one parameterization to another. The parameterization-dependence of shell effects becomes obvious when one directly compares the entire barriers, as can be seen from Fig. 6. Even the overall shape of the barrier is not the same for all parameterizations. The deformation of the fission isomer varies, as does the deformation of the configuration corresponding to the top of the barriers. Even more intriguingly, for some parameterizations there appears a third minimum at large asymmetric shapes, or the topography around the spherical point is quite different. For the purpose of the present paper, it is not important to disentangle the origin of these variations, which are for example related to the strength of spin-orbit and tensor terms as will be discussed in a forthcoming article. What is relevant is that there obviously is a large variation of the deformation dependence of the shell effects among the parameterizations studied here. This, in turn, indicates that the adjustment of the excitation energy of the fission isomer or the fission barrier heights cannot be easily replaced by an adjustment of a universal empirical value of the surface energy coefficient.  [69] and SLy5* [70] parameterizations are also given for comparison. The first block (columns 2-8) shows the standard bulk properties for infinite symmetric matter, i.e. saturation density ρsat in fm −3 , energy per particle E/A in MeV, incompressibility K∞ in MeV, symmetry energy coefficient J and its slope L in MeV, and isoscalar effective mass m * 0 /m and enhancement factor of the TRK sum rule κv. The second block (semi-bulk; columns 9-11) shows surface energy coefficients a surf in MeV as obtained from semi-infinite nuclear matter calculations within the MTF, HF and ETF schemes (see text This is complicated further by the apparent impossibility of determining a precise value for a surf in a model-independent way. As indicated by the analyses of Refs. [10,33], a reliable extraction of a surf and other liquid-drop parameters from microscopically calculated binding energies of finite nuclei requires to go to systems with A ≈ 10 5 nucleons. The same can be expected to hold for its reliable extraction from the binding energies of real nuclei, but the systems needed do not exist in nature. Second, the preferred value for a surf will also depend on choices made for quantum corrections to the binding energy in a given fit. For example, the size of the rotational correction increases rapidly with deformation [1,21,99], such that parameterizations that are supposed to be used with it require a larger value for a surf than parameterizations supposed to be used without it. Also, for parameterizations that are supposed to be used in beyond-mean-field models, it is the energy difference between collective states that should be compared with data [99], not the difference of minima of the energy surface.

VI. A FIT PROTOCOL INCLUDING a surf
In one way or the other, information about fission barriers has already been sometimes used to constrain parameterizations of effective interactions of Skyrme [7,18,19,21] type.
Using the fission barriers themselves for the parameter adjustment, however, is quite costly. Assuming that shell effects are principally fixed by the other ingredients of the fit protocol, our results suggest a way to fit the information contained in fission barriers in a more economical iterative process: A parametrization is first adjusted to reproduce an initial guess for the value of a surf , or a series of values of a surf that cover a reasonable range. After convergence of the fit, one calculates then a couple of well-selected fission barriers with the preliminary parameterizations, checks the deviation from data, estimates how a surf should be changed and runs a new fit (or series of fits) with the improved estimate for a surf . The process can be repeated until the desired quality for fission barriers is reached. Figure 5 suggests that the correlation between the surface energy energy coefficient and the characteristic features of the energy surface is basically the same within each method to calculate a surf , in spite of the HF, ETF and MTF calculations providing different values for a given parameterization. This means that any of these methods can then be used during the fit as long as the value used for a surf are tuned accordingly to the formalism used.
Because of its computational simplicity, we will use the MTF value for a surf to construct a series of fits with a surface energy coefficient systematically varied in the range between 18.0 and 19.4 MeV in steps of 0.2 MeV, while everything else in the fit protocol is kept unchanged. As a starting point, we chose the protocol used to adjust SLy5* [70], which differs from the fit protocol of the original SLy5 parameterization [69] mainly by an additional constraint that avoids the appearance of unphysical finite-size instabilities in the spin channels [70]. In order to have a series of parameterizations with similar bulk properties, we have added an additional constraint on the value of the L coefficient of nuclear matter, which otherwise is only scarcely constrained by finite nuclei, such that it remains in the interval of (50 ± 2) MeV.  Table V. When comparing any two parameterizations, the difference between the values of the surface energy coefficients obtained with different methods is almost constant, see Fig. 7. Figure 8 displays the deformation energy along the static fission path of 240 Pu obtained with these parameterizations. As they are all fitted within the same protocol but for the value a surf is fixed at, the evolution of the curves with a surf is now much more regular. However, there are still indications that also shell effects are slowly varying in response to the change of a surf in the fit pro-  tocol: at β 2 values around 1.7, a shallow third minimum develops with increasing value of a surf .
Still, for the excitation energies of the fission isomer and the inner and outer saddles, the correlation with a surf is now almost perfectly linear as demonstrated by Fig. 9. Still, none of the parameterizations does simultaneously describe all of the three properties.
In the present study, we have focused on the isoscalar surface energy coefficient a surf . All candidates for heavy nuclei whose deformation energy can be used to constrain parameter fits, however, have a neutron excess. For example, the nucleus 240 Pu used as a benchmark above has an asymmetry of I ≡ (N − Z)/(N + Z) = 0.2. The leptodermous expansion of the nuclear binding energy suggests that there is a correction to the surface energy that depends explicitly on the nucleus' asymmetry I and that in the liquid-drop model is parametrized through the surface symmetry energy coefficient a ssym [22,33,62,63].
If a single nucleus is used to constrain the surface energy, there is the danger that a ssym accidentally takes a wrong value and thereby introduces an erroneous isospin dependence of fission barriers. However, as demonstrated in a detailed analysis of this quantity's role for the sys-tematics of the excitation energy of superdeformed states of heavy nuclei [22], the value of the a ssym of existing parameterizations of the Skyrme EDF is strongly correlated with the volume energy coefficient and therefore fixed by masses along the valley of stability, such that it cannot vary freely over a wide interval. This gives us confidence that a ssym takes a reasonable and consistent value for all SLy5sX fits. Indeed, for all of the SLy5sX parameterizations we find values of a ssym that are close to −49 MeV. Also, calculations of the fission barriers of other heavy nuclei that will be reported elsewhere indicate that the overall trend of the predicted fission barriers does not change with asymmetry. For example, when calculated with SLy5s1, the fission barrier height of the much less asymmetric nucleus 180 Hg with I = 0.11 is described as well as is the one of 240 Pu.

VII. DISCUSSION AND CONCLUSIONS
To summarize our main findings concerning the various methods to calculate a surf : (i) HF, ETF and MTF calculations for semi-infinite nuclear matter provide very consistent, but not identical, values for a surf . The MTF method gives always slightly larger values than the HF approach, whereas the ETF method always gives slightly smaller values than HF. Differences between the values extracted from semi-classical MTF and ETF calculations on the one hand and the quantal HF approach on the other hand are typically on the order of 500 keV out of about 18 MeV, but in exceptional cases might be as large as 1 MeV.
(ii) The values of a surf as extracted from ETF calculations of semi-infinite nuclear matter and from the systematics of the ETF binding energies of very large artificial spherical nuclei are very close and differ rarely by more than 10 keV. To reach this level of agreement, an A 0 term has to be included in the leptodermous expansion of the binding energy of finite nuclei, and the contribution of the center-of-mass correction omitted.
(iii) For the purpose of calculating a surf , the MTF and ETF approximations are fairly robust. The deviation from HF values is reasonably parameterization independent, although there are differences in detail that are correlated to the presence and strength of the tensor terms, the incompressibility K ∞ , and the isoscalar effective mass m * 0 /m. During a fit, the values for K ∞ , m * 0 /m and the coupling constants of the tensor terms will rarely vary over a large interval, such that the analytical MTF value for a surf as obtained from from Eq. (22) can be safely used in a fit protocol to determine the coupling constants of Skyrme EDFs. In any event, we do not aim at the reproduction of a universal empirical value for a surf , which will be difficult to extract in a modelindependent way from data anyway, but instead provide a simple and efficient control over the surface properties within a given framework.
Concerning the correlation between characteristic energies of the fission barrier of 240 Pu and the values for a surf we find that (i) The simultaneous description of the fission barrier heights and the excitation energy of the fission isomer is not trivial as they are also strongly sensitive to shell effects which are at the very origin of the complicated topography of the deformation energy surface.
(ii) While the characteristic energies of fission barriers are clearly and unambiguously correlated to the surface energy coefficient in the expected manner, for existing Skyrme parameterizations there is a large scatter when plotting them as a function of the surface energy coefficient. This is a consequence of the shell effects being unsystematically different when parameterizations are constructed with different fit protocols.
To eliminate the protocol-dependence of the analysis of a surf and its correlation with fission barriers, and as a proof of principle for a parameter fit that includes the MTF value for a surf in its protocol, we have constructed a series of parameterizations of the standard Skyrme functional called SLy5s1, SLy5s2, . . . , SLy5s8, with systematically varied surface energy coefficient.
(i) The resulting parameterizations exhibit almost linear correlations between the surface energy coefficient on the one hand and each of the characteristic energies of the energy surface of 240 Pu on the other hand, providing the proof-of-principle for replacing the adjustment of fission barriers with a suitably chosen value of a surf as obtained from semi-infinite nuclear matter calculations with any of the schemes considered here.
(ii) None of the SLy5sX parameterizations, however, reproduces simultaneously the empirical data for the excitation energy of the fission isomer as well as the heights of the inner and outer barriers. This points to imperfections of the deformationdependence of shell effects provided by these parameterizations.
(iii) The best reproduction of the energy surface of 240 Pu is provided by SLy5s1, the parameterization with the lowest value of a surf considered in the series of fits.
Adjusting parameters to reproduce a value for a surf is computationally much easier than reproducing fission barriers of heavy nuclei, not only with respect to CPU time, but, even more importantly, also in terms of stability and controllability of the calculations during the early stages of a fit when the parameters are still far from their physical values such that the topography of the energy surfaces might be very different from the physical one. In practice, however, the appropriate value of a surf that provides the best description of the targeted deformation energies will have to be determined iteratively, alternating between complete parameter fits for some guess(es) for a surf on the one hand, and the calculation of deformation energy surfaces with the such determined parameterizations on the other hand. The latter then provide a refined guess for a surf , if necessary. In the present paper, we limited ourselves to the first of such cycles. As many of the widely-used Skyrme parameterizations provide an unsatisfactory description of the systematics of fission barriers and excitation energies of very deformed states [1,22,101,111], it is highly desirable to better constrain the nuclear surface energy in future fits. Besides the obvious importance for the study and understanding of the fission process itself, fission barrier heights are also a determining factor for the stability of superheavy nuclei [101,111] and the dynamics of the astrophysical r-process of nucleosynthesis [21,49,111].
In conclusion, the use of the MTF pocket formula for a surf provides a rapid and robust expression to control the surface properties of nuclear energy density functionals during the adjustment of their parameter. In this paper, we have focused on the (isoscalar) surface energy coefficient. Within the MTF approach, the integral that appears in the calculation of the (isovector) surface symmetry energy coefficient from semi-infinite matter cannot be solved analytically. However, making further approximations whose consequences remain to be analyzed, one can arrive at an analytical expression for this quantity as well [112], which might offer a route to its efficient fine-tuning.
The generalization of the MTF pocket formula to Skyrme-type functionals with derivative terms of order 4 or 6 [13][14][15] is feasible, but will require an arduous extension of the formalism to higher-order terms.