{'fraction_non_alphanumeric': 0.0336177820373408, 'fraction_numerical': 0.016772029415559282, 'mean_word_length': 4.733093829247675, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nuclear energy density optimization: Shell structure

Nuclear density functional theory is the only microscopical theory that can be applied throughout the entire nuclear landscape. Its key ingredient is the energy density functional. In this work, we propose a new parameterization UNEDF2 of the Skyrme energy density functional. The functional optimization is carried out using the POUNDerS optimization algorithm within the framework of the Skyrme Hartree-Fock-Bogoliubov theory. Compared to the previous parameterization UNEDF1, restrictions on the tensor term of the energy density have been lifted, yielding a very general form of the energy density functional up to second order in derivatives of the one-body density matrix. In order to impose constraints on all the parameters of the functional, selected data on single-particle splittings in spherical doubly-magic nuclei have been included into the experimental dataset. The agreement with both bulk and spectroscopic nuclear properties achieved by the resulting UNEDF2 parameterization is comparable with UNEDF1. While there is a small improvement on single-particle spectra and binding energies of closed shell nuclei, the reproduction of fission barriers and fission isomer excitation energies has degraded. As compared to previous UNEDF parameterizations, the parameter confidence interval for UNEDF2 is narrower. In particular, our results overlap well with those obtained in previous systematic studies of the spin-orbit and tensor terms. UNEDF2 can be viewed as an all-around Skyrme EDF that performs reasonably well for both global nuclear properties and shell structure. However, after adding new data aiming to better constrain the nuclear functional, its quality has improved only marginally. These results suggest that the standard Skyrme energy density has reached its limits and significant changes to the form of the functional are needed.


I. INTRODUCTION
An important goal in research in low-energy nuclear physics is to develop an universal nuclear energy density functional (EDF) that can be used to explain and predict static and dynamic properties of atomic nuclei within the framework of nuclear density functional theory (DFT). Building such a functional has been one of the primary drivers behind the formation of the former UNEDF SciDAC-2 collaboration [1,2]; its current successor, the NUCLEI SciDAC-3 collaboration [3]; and the FIDIPRO collaboration [4].
Most of the nuclear EDFs used in self-consistent meanfield calculations have been derived from phenomenological effective interactions, or pseudopotentials [5][6][7]. A recent promising development is to use results from effective field theory combined with density matrix expansion techniques to construct a realistic EDF based on chiral interactions [8][9][10][11][12]. In a parallel effort, methodologies have been developed to validate nuclear EDFs by optimizing their low-energy coupling constants to experimental data on finite nuclei and pseudo-data on nuclear matter and other relevant systems [2,[13][14][15][16]. One of the main challenges of EDF optimization is to find the most relevant fit observables that can tightly constrain the parameter space of the model. This requires a careful analysis of the contribution of each term of the functional to low-energy nuclear properties. In finite nuclei, the most common observables used in EDF fits are binding energies and their differences, charge radii, surface thickness, and energies of giant resonances [7,17]. Given a mathematical form of the EDF, its predictive power ultimately depends on the choice of the data used in the optimization. In particular, applications of DFT to nuclear spectroscopy are very sensitive to the details of shell structure; this requires a careful choice of fit observables.
Shell structure is the fundamental property of the atomic nucleus [18]. In an independent-particle picture, arXiv:1312.1746v3 [nucl-th] 28 Apr 2014 shell structure can be associated with the single-particle (s.p.) spectra of the mean-field potential [19,20]. Reproducing the correct ordering and distribution of s.p. levels is, therefore, an essential requirement for nuclear structure theories, but it has to be approached with caution [19,21,22] since s.p. motion is significantly modified by correlations [22][23][24]. In the context of nuclear DFT, many commonly used EDFs have been optimized by explicitly using some experimental input pertaining to the s.p. level structure in doubly-magic nuclei [7].
The s.p. shell structure is very sensitive to the details of the effective interaction or the energy density and is the result of a subtle interplay between the gradient terms and effective mass, spin-orbit, and tensor terms [25,26]. Suggestions to study tensor interactions within the selfconsistent mean-field approach were made already in the seventies [27] but the limited experimental data available did not provide sufficient sensitivity to adjust the related coupling constants. In recent years, the role of tensor coupling constants, in Skyrme EDFs in particular, has been thoroughly investigated [25,26,[28][29][30][31][32][33][34][35][36][37][38]. An important conclusion from several of those papers is that the inclusion of tensor terms should not be done perturbatively but should instead involve the complete EDF reoptimization at the deformed Hartree-Fock-Bogoliubov (HFB) level. This implies that constraints on the tensor terms must be included in the pool of fit observables.
In our previous works on Skyrme EDF optimization [14,15], tensor terms were disregarded because our dataset did not contain any information specifically constraining shell structure. This limitation is lifted in this article, which should be viewed as the continuation of our work on energy density parameter optimization. In Ref. [14], we presented the main strategy underlying our optimization protocol and developed the unedf0 EDF parameterization by using experimental input on a selected set of nuclear masses, charge radii, and odd-even mass differences. In the same paper, we performed one of the first sensitivity analyses of Skyrme parameterizations to obtain the correlations and standard deviations for the parameters. In Ref. [15], we modified the form of the functional by removing the center of mass correction, which allows for straightforward time-dependent Hartree-Fock (HF) and HFB applications. To constrain deformation properties, we also extended our dataset to include information on fission isomer excitation energies. In this way, nuclear deformation properties produced by the resulting unedf1 functional have greatly improved, in particular in the context of nuclear fission.
The goal of this study is to include the tensor coupling constants in the set of optimized parameters. The resulting energy density is a general functional of the one-body density matrix up to second order in derivatives. Because shell structure is very sensitive to the tensor terms of the functional, we extend our experimental dataset by adding a set of s.p. energy splittings in doubly-magic nuclei. The optimization of the functional within this extended dataset yields the unedf2 parameterization of the Skyrme energy density. In the spirit of our previous work, we carry out a full sensitivity analysis of unedf2, which is essential to assessing the predictive power of the theory [39][40][41][42][43][44][45][46]. Since the previous work of Ref. [26] based on the linear regression methodology demonstrated that the current standard form of the Skyrme EDF cannot ensure a spectroscopic-quality description of s.p. energies, unedf2 is certainly not the universal nuclear EDF. It can be viewed, however, as the best all-around Skyrme EDF that performs reasonably well for both global nuclear properties and shell structure. For this reason, we consider unedf2 as the end of the standard Skyrme EDF journey.
This paper is organized as follows. In Sec. II, we briefly review the theoretical framework and the notations. Section III describes the optimization method employed, experimental data used in the fit, and presents the unedf2 parameterization together with its sensitivity analysis. Global nuclear properties computed with unedf2 are reviewed in Sec. IV. Finally, Sec. V contains conclusions and perspectives for future work.

II. THEORETICAL FRAMEWORK
In the nuclear DFT, the total energy E is a functional of the one-body density matrix ρ and pairing densityρ and can be cast into the generic form where E Kin (r) is the kinetic energy; χ t (r) is the isoscalar (t = 0) and isovector (t = 1) particle-hole Skyrme energy density;χ(r) is the pairing energy density; and E Coul (r) is the Coulomb term. The particle-hole part of the Skyrme energy density reads where each term is multiplied by a coupling constant C uu t represented by a real number. The coupling constant C ρρ t is the only exception, as it has the traditional densitydependence The definitions of the various densities ρ, τ , and J µν (J is the vector part of J µν ) can be found in Ref. [7,30,[47][48][49].
The coupling constants C ρρ t and C ρτ t are related to the volume part of the energy density and can be expressed as a function of the parameters of infinite nuclear matter [14].
The term µν J µν,t J µν,t in Eq. (2) represents a tensor energy density. The approximation that was made in the unedf0 and unedf1 optimizations was to set C JJ t = 0 for both t = 0 and t = 1. In the present work, this constraint has been removed, and these two tensor coupling constants are taken as free parameters. Note that we do not allow independent variations of the pseudoscalar, vector, and peudotensor components of J µν,t : each of these components is multiplied by the same coupling constant, see [48,49].
The pairing term is derived from the mixed pairing force of Ref. [50], leading to the energy densitỹ where V q 0 (q = n, p) is the pairing strength. In this work as before, we take ρ c =0.16 fm −3 . We have used different pairing strengths for neutrons and protons [51]. Owing to the zero range of our effective pairing force, we have used a pairing cut-off E cut = 60 MeV to truncate the quasi-particle space. To prevent the collapse of pairing correlations near closed shells, we have also used the variant of the Lipkin-Nogami (LN) method as in Ref. [52].
As with unedf1, we disregard the center-of-mass correction; for motivation see Ref. [15]. The Coulomb term contains direct and exchange contributions. The direct part is computed from the proton density distribution assuming a point-proton charge, and the exchange part is treated at the Slater approximation.

III. OPTIMIZATION OF ENERGY DENSITY AND SENSITIVITY ANALYSIS
In this section, we present our optimization protocol and analyze the features of the resulting unedf2 parameterization. In Sec. III A, we review the experimental data used in the fit. In Sec. III B, we describe how information on shell structure was incorporated. The optimization procedure itself is summarized in Sec. III C. Section III D compares unedf2 to previous unedf parameterizations, and Sec. III E contains the results of the sensitivity analysis.

A. Experimental Dataset
In order to determine the tensor coupling constants C JJ t , one must properly select the experimental fit observables to effectively constrain their values. To this end, we have extended the previous dataset used in the unedf1 optimization by including an additional nine single-particle level splittings, five new data points for odd-even staggering (OES), and one additional binding energy.
We show in Table I the empirical values of singleparticle splittings in several doubly-magic nuclei. All values are taken from the empirical s.p. energies listed in Ref. [53] except for the proton f 5/2 − f 7/2 splitting in  48 Ca, which is taken from Ref. [54]. The rationale to use s.p. splittings instead of the absolute energy of levels is to remove some of the systematic errors induced by the use of a truncated harmonic oscillator (HO) basis. We set the weight of the single-particle data points in the χ 2 function to w = 1.2 MeV. This choice was motivated based on the singular value decomposition (SVD) analysis performed in Ref. [26], which showed that Skyrme EDFs can reproduce empirical s.p. levels at this precision level. We recall that the weight can be viewed as a coarse estimate of the theoretical error on a given observable. The calculation of s.p. splittings in 132 Sn requires the ground-state energy of 132 Sn, see Sec. III B. For this reason, we added the binding energy of this nucleus to the dataset. As in Ref. [14], experimental information has been taken from the 2003 mass evaluation, and the nuclear binding energy was obtained after taking into account the electronic correction, yielding the value B( 132 Sn) = −1102.686066 MeV. For this additional datum, we took the same weight w = 2 MeV as for other binding energies. In addition to the single-particle splittings and the binding energy of 132 Sn, we have added five new OES data points, which are listed in Table II. This was motivated by the observation that pairing properties of actinide nuclei and neutron-rich tin isotopes are poorly reproduced by unedf1, suggesting that the weight of pairing-related data in the objective function should be increased. We recall that the experimental OES that we use is defined as the average of two odd-even (protons) and even-odd (neutrons) ∆ To summarize, the unedf2 optimization dataset contains 47 deformed binding energies, 29 spherical binding energies, 28 proton point radii, 13 OES values, 4 fission isomer excitation energies, and 9 single-particle level splittings. The changes with respect to the unedf1 optimization procedure are as follows: • Tensor coupling constants C JJ 0 and C JJ 1 are optimized; • The binding energy of 132 Sn is added to dataset; • 9 new single-particle splittings are included in dataset; • 5 new OES data points are added to dataset; • The weight of all OES data points is increased to 0.1 MeV.

B. Computation of Single-Particle Levels
Most of the optimizations of the Skyrme EDF that included information on s.p. splittings were conducted at the spherical HF level (see, for example, Refs. [13,55,56]). Within this approximation, the many-body wavefunction reduces to a single Slater determinant, and the theoretical single-particle levels are taken as the eigenvalues of the HF Hamiltonian ε j following Koopmans' theorem [26,57]. In the framework of the HFB theory with approximate particle number projection, where the basic degrees of freedom are not particles but quasi-particles, separation energies provide a more convenient quantity to relate to effective single-particle energies.
In our optimization procedure, theoretical singleparticle splittings were thus obtained by applying the blocking HFB approach with the approximate LN correction. We employ the equal-filling approximation to blocking, since it yields an excellent estimate of the full symmetry-breaking blocking results [58]. Since the shape polarization induced by the blocking prescription spontaneously breaks the spherical symmetry of the odd-A nucleus, the energy degeneracy of a blocked spherical quasi-particle orbital with angular momentum j is lifted [33,58,59]. In the equal filling approximation, however, time-reversal symmetry is conserved, and states with the angular momentum projection +Ω and -Ω are degenerate.
This fragmentation of any given j-shell in nonspherical blocking calculations poses a practical difficulty. Indeed, one should in principle compare the energy of all obtained blocking configurations with different |Ω| values, and pick the lowest one to compare with experiment. The difficulty with such a strategy is that a configuration with |Ω| < j can originate from a j-shell that is different from the one under consideration. To avoid such a situation, we have chosen to block the single state with the maximum projection Ω = +j. The associated systematic error does not exceed 100 keV [15].
Empirical s.p. energies are usually extracted from the centroids of (often broad) strength functions of pickup/stripping reactions. In our approach, we choose to relate these empirical levels to one-particle separation energies computed at the HFB+LN approximation. This choice allows us to remain consistent throughout and calculate all observables at the same approximation level (HFB+LN). In addition, we find that empirical energy splittings extracted from the s.p. energies of the eveneven nucleus or directly from the separation energy of the odd nucleus in a given J π configuration differs typically by at most a few hundreds keV, with one notable exception. This should be compared with the >1.2 MeV accuracy of Skyrme functionals for s.p. data, and suggests that the determination of the empirical value should not have a large impact on the optimization. The exception for the neutron f 5/2 state in 40 Ca is due to the strong fragmentation of the f 5/2 strength among multiple states, resulting in a rather broad centroid [60]. Similarly, the proton f 5/2 state in 48 Ca is fragmented among multiple states [61]. In the Supplemental Material at [62], we provide both sets of experimental s.p. splittings for convenience.
Our procedure to generate theoretical s.p. splittings follows that of Ref. [63] and can be summarized as follows. We begin by computing a reference spectrum in the doubly magic nucleus of interest. Next, we use this quasiparticle spectrum to identify the blocking configuration with Ω = +j and perform the blocking calculation in the system with ±1 particles. The effective particle and hole s.p. energies are respectively defined as where A refers to the particle number of the reference (doubly-magic) nucleus and E bl is the energy of the blocked configuration in the neighboring odd nucleus. The labels "hole" and "particle" refer to whether the corresponding s.p. levels would be, respectively, fully occupied or empty in the corresponding HF calculation of the doubly-magic nucleus. In the case where the s.p. levels involved in the s.p. splitting are both either above or below the Fermi surface, the contribution of the even-even binding energy cancels out, and the s.p. splitting reduces to the difference of total binding energies of blocked configurations.

C. Optimization
All HFB calculations in the unedf2 optimization were performed with the DFT solver hfbtho [64]. The code solves the HFB equations in an axially symmetric deformed HO basis. In our initial work on unedf0, we used a spherical HO basis with 20 shells and assumed the HFB solution to be reflection symmetric. Adding experimental data on fission isomer excitation energies in the unedf1 optimization required computing the energy of super-deformed (SD) configurations. In order to mitigate truncation errors, the SD states were calculated with a deformed, or stretched, HO basis with the axial quadrupole deformation parameter β = 0.4. In this work, we have maintained the same setup as for unedf1: the spherical basis is used for all ground-state configurations, and the stretched basis with β = 0.4 is used for SD states. In all cases, the spherical frequency ω 0 of the HO basis is set at ω 0 = 41/A 1/3 [64].
The objective function in our optimization is where D T is the number of different data types; for unedf2 D T = 5. The total number of data points is , and the number of parameters to be fitted is n x (here n x = 14). The calculated value of the j th observable of type i is s i,j (x), while the corresponding experimental value is denoted by d i,j . Each data type has a weight w i . As in the cases of un-edf0 and unedf1, the isovector effective mass 1/M * v is kept constant during the optimization, since we find that it cannot be constrained reliably with the current data; we retain the SLy4 value for 1/M * v for historical reasons [56]. In the Supplemental Material [62], we provide the full list of experimental data points that were used in all three optimizations: unedf0, unedf1, and unedf2.
As in our previous work, the parameters of the functional are not allowed to attain unphysical values: we impose bounds on the range of variation of each parameter. Bounds for the parameters common to both unedf1 and unedf2 were assumed to be the same. We did not set any bounds on the tensor coupling constants C JJ t . During the optimization, two of the parameters, E/A and L, ran to their boundary and were fixed to those values. They were subsequently excluded from the sensitivity analysis. The optimization was carried out with the same pounders algorithm that was used for other unedf parameterizations, see Ref. [14] for details.

D. The UNEDF2 parameterization
The optimized parameter set of the EDF unedf2 is listed in Table III along with the standard deviation of each parameter and the 95% confidence intervals. To facilitate legibility, we display only the first few significant digits of each parameter. In the Supplemental Material [62] we provide the parameter values of all three parameterizations up to machine precision in two different representations: the hybrid nuclear matter/coupling constants representation and the full coupling constant representation. We first note that for unedf2, the nuclear incompressibility parameter K, while in the top range of acceptable values, is now constrained by the data, whereas the slope of the symmetry energy L is not. As expected, both the neutron and pairing strengths are also a little larger, a direct consequence of adding more OES points into the dataset. Table IV lists all three unedf parameterizations produced so far, and compares them to the SLy4 parametrization, which was the starting point for un-edf0. Interestingly, the unedf2 and unedf1 parameterizations are quite similar overall. This result is a little surprising: one may have expected that relaxing the constraints on the tensor coupling constants would lead to a significant rearrangement of all other coupling constants, in particular the spin-orbit coupling constants. Indeed, it was shown in Ref. [30] that there is a strong anti-correlation between the isoscalar spin-orbit and tensor coupling constants. This relationship is confirmed in our optimization through a large correlation coefficient of −0.88 between C ρ∇J 0 and C JJ 0 . In fact, the values of C ρ∇J 0 and C JJ 0 are consistent with the empirical C ρ∇J 0 (C JJ 0 ) dependence reported in Ref. [30]. Yet, in spite of this very strong correlation, the value of C ρ∇J 0 changes only by 13% between unedf1 and unedf2.
Looking more closely at the values of all spin-orbit and tensor coupling constants of the unedf2 parameterization, we find that they are compatible with the results of Ref. [33]. In particular, the value of C ρ∇J 0 is close to the "universal" value of −60 MeV obtained from a refit of three different EDFs to the spin-orbit splittings in 40,48 Ca and 56 Ni. Our optimized tensor coupling constants, as well as their isoscalar-isovector trend, are also in the same ballpark as those partial refits of Ref. [33]. They are, however, just beyond the space of the TIJ family of parameterizations considered in Ref. [30]. Recent work also suggests that only the region where C JJ 0 + C JJ 1 < 0 and C JJ 0 − C JJ 1 > 0 should be physical [38]. This is the case for the unedf2 functional. Since the unedf2 fit is an optimization carried out by considering a broad range of nuclear properties (with five different types of experimental data), it is encouraging that our results overlap well with those obtained in systematic studies of spin-orbit and tensor terms. We show in Table V the unedf2 coupling constants in natural units [65,66]. According to the hypothesis of naturalness, the magnitude of (the absolute value of) coupling constants should be of order unity, when scaled into unitless quantities. The scale Λ used to perform the transformation to natural units was taken as Λ = 687 MeV, which was found in Ref. [66] to be valid for Skyrme EDFs. As seen in Table V, nearly all the unedf2 coupling constants fall in the interval [1/3, 3] which is compatible with the hypothesis of naturalness [66]. The one notable exception is C ρτ 0 , which is unnaturally small; C ρρ 1D and C ρτ 1 are also at the limits of the allowed interval.

E. Sensitivity Analysis
The standard deviations σ of the unedf2 parameterization are listed in Table III, together with the 95% confidence intervals. We recall that the standard deviations (and also correlations) are calculated only among those parameters that do not run into their imposed boundaries. Compared with the previous parameterizations unedf0 and unedf1, the standard deviations are overall smaller, reflecting improved constraints on the coupling constants. For example, the standard deviation of the symmetry energy went down from 3.05 MeV for unedf0 to 0.60 MeV for unedf1 to only 0.32 MeV for unedf2. Similarly, the isoscalar effective mass, which could not be constrained in unedf0, had a standard deviation of 0.12 for unedf1, which was further reduced to 0.05 for unedf2. This improvement on constraining all coupling constants of the functional, while not perfect, is a confirmation of the validity of our strategy.
Table VI displays the correlation matrix among the coupling constants of the functional. As with unedf1, there exists a strong correlation between the pairing strength parameters and both the isoscalar effective mass and the isoscalar surface coefficient C ρ∆ρ 0 , although for the latter the correlation is less pronounced than for unedf1. These correlations reflect a strong interplay between the level density near the Fermi level and the magnitude of pairing correlations. Another strong (anti-)correlation can be observed between ρ c and K. This is reminiscent of our first parameterization, dubbed unedfnb in Ref. [14], where no bounds had been imposed on the coupling constants. Further, as discussed in Sec. III D, the isoscalar spin-orbit and tensor coupling constants are strongly anti-correlated and seem to follow the trend predicted in Ref. [30].  The overall impact of each data type on the unedf2 parameterization can be assessed by studying the sensitivity matrix S: where J(x) is the Jacobian matrix calculated with the Following Refs. [14,15], we have calculated the partial sums of the absolute values in S(x) for each data type, normalized with respect to the number of data points. The results are presented in Fig. 1, with each bar normalized to 100%. A number of observations made for the unedf0 or unedf1 functionals still apply, such as the strong sensitivity of pairing strengths to OES data or the large impact of fission isomer excitation energies on the determination of a sym . Overall, s.p. splittings, fission isomer excitation energies, and OES data seem to be the main drivers of the parameterization, while the relative role of masses is reduced. Looking closely at the coupling constants that are relatively well constrained, one may identify two trends: (i) bulk coupling constants (i.e., ρ, K, and a sym ) are not really impacted by the OES data; (ii) surface coupling constants (involving gradient terms) are more sensitive to OES data, fission isomer excitation energies, and s.p. splittings. The three isovector surface coupling constants (C ρ∆ρ 1 , C ρ∇J 1 , C JJ 1 ) behave differently but are less constrained by the data, as shown by their large standard deviations reported in Table III. A complementary way to examine our optimization dataset is to analyze the impact of individual data points on the optimized solution. This is plotted in Fig 2. Here, the amount of variation for the optimal solution is presented when each data point d i,j is shifted by an amount of 0.1 w i . The standard deviations σ k for parameters x k are listed in Table  III. As for unedf0 and unedf1, the overall changes in x are of the same order of magnitude and are very small, δx/σ ≈ 0.01. This indicates that the set of fit observables in unedf2 has been chosen consistently. The new s.p. data points seem to have a relatively large impact on the parameterization, especially the s.p. splittings in 40 Ca.

IV. PROPERTIES OF UNEDF2 PARAMETERIZATION
In this section, we review various properties of the un-edf2 parameterization. In Sec. IV A, we apply the linear response theory to test the functional against the presence of finite-size instabilities. Section IV B discusses correlations among various observables. Predictions of unedf2 for shell structure in doubly-magic nuclei are presented in Sec. IV C, and global binding energy and deformation trends (in particular in the context of nuclear fission) are analyzed in Secs. IV D and IV E, respectively. Section IV F contains the discussion of unedf2 predictions for neutron droplets in external traps.

A. Linear Response and Instabilities
The linear response formalism in nuclear physics has been developed mainly in the framework of the Random Phase Approximation (RPA) based on the use of an effective interaction. In Refs. [67] this formalism was generalized to determine the response function in both symmetric nuclear matter and pure neutron matter for the case of a general Skyrme EDF as given in Ref. [49]. This recent development is needed here because of the nonstandard spin-orbit term of the unedf parametrizations. The response functions χ (S,M,T ) (ω, q) of interest to the present work are defined as the response of the infinite medium to external probes of the typê where S (M ) is the spin (its projection along the z-axis), T is the isospin, and Θ

(S,M,T ) j
is an operator acting on spin, isospin, or both. As recently shown in Ref. [68], the response function of neutron matter can provide information about instabilities in finite nuclei [69]. More precisely, it was shown that whenever a pole appears in the response function close enough to the saturation density of the system ρ c , the finite nucleus undergoes an instability in the corresponding channel (scalar/vector, isoscalar/isovector), see examples in Ref. [58].
Several quantitative criteria to estimate the likelihood of finite-size instabilities for a given EDF have been recently proposed [70,71]. Because of shell fluctuations, the nucleus can explore regions of densities slightly larger than the saturation density. In Ref. [71], the following conservative criterion was established: whenever the response function has a pole at a density ρ ≈ 1.4ρ c , there is a risk of instability in calculations for finite nuclei. More complex criteria were proposed in [67]. The poles of the response function are in practice determined by solving the equation where the expressions of χ (S,M,T ) can be found in Ref. [67]. Figure 3 shows the position of the lowest poles of the response function as a function of the transferred momentum q in symmetric nuclear matter for unedf0, un-edf1, unedf2, and SkP [72]. Since the unedf functionals have been developed to be used only in the time-even channel, we will limit our analysis to S = 0. We first remark that in all cases we have an instability in the region at low density and low momentum in the channel S = 0, M = 0, T = 0. This is the well-known spinodal instability, which is physical. We then observe that all unedf functionals satisfy the stability criterion given in Ref. [71]. On the other hand, SkP is unstable in the scalar/isovector channel in the region of densities around ρ c ; as a result it was shown in Ref. [69] that calculations of finite nuclei with this functional are unstable. Although a quantitative criterion is not yet available for pure neutron matter, it is interesting to study the instabilities of such a system. Figure 4 shows the position of the lowest poles for pure neutron matter. While the unedf0, unedf1, and (to a lesser extent) SkP EDFs do not exhibit poles near the saturation density, the situation is different for unedf2. This suggests that S = 0 instabilities could manifest in neutron-rich nuclei or in trapped neutron droplets. As we will see in Sec. IV F below, such instabilities are indeed present in heavy neutron droplets calculated with unedf2. This result may be a consequence of the large negative values of the tensor coupling constants [30].

B. Correlations with Other Observables
The sensitivity analysis presented in Sec. III E aimed at quantifying the behavior of the χ 2 landscape at the minimum for the set of observables used in the fit. Complementary information can be obtained from an analysis of the correlations between observables not included in the fit [39,42,43]. These correlations can be extracted from an estimate of confidence regions near the minimum. In this section, we define the confidence region based on a criterion for the value of the objective function [73]. It is a slightly different prescription from the procedure that we use to define the confidence intervals; see Sec. III.B.1. of Ref. [14] for details. Asymptotically (at the limit of n d → +∞), both prescriptions are in fact equivalent (see discussion in Sec. 3.3.1 of Ref. [74]).   Here we construct an approximate confidence region using the following approach. At the minimum x fin , the quantity χ 0 = χ 2 (x fin ) characterizes the best-fit parameterization. The parameters x in a neighborhood V of the minimum can still provide a reasonable description of nuclear properties. The parameter space V is thus referred to as a "reasonable" domain. Each observable A that can be computed in the EDF theory is also a function of the Skyrme parameters, A = A(x). Varying x in the vicinity of the optimal set will lead to fluctuations in the values of A with respect to its value at the minimum, A 0 = A(x fin ). The uncertainty of the prediction is characterized by the variance Var(A) = σ 2 A = (A − A 0 ) 2 , where the average value is computed from This simple estimate of uncertainties provides valuable information on the predictive power of the model. Fur-ther information can be obtained from the correlation coefficient c AB between two observables A and B defined from the covariance matrix as

Var(A)Var(B)
.  Table VII. Table VII shows the values and the uncertainties of a large set of observables. The vicinity V was defined by the level set V = x : χ 2 (x) − χ 0 ≤ n d − n x . As noted earlier, this construction of V is different from the one used to define the 95% confidence interval, and hence, the standard deviations reported in Table VII are slightly  different from those of Table III. As discussed later in Sec. III E, the set unedf2 delivers rather small uncertainties for all observables shown. Figure 5 shows the correlation matrix c AB (12) between various pairs of observables computed for 208 Pb. To illustrate the impact of the optimization protocol, we compare unedf2 with the SV-min parameterization [13]. In the case of SV-min (upper panel), we can see four blocks of highly correlated observables [44,46,75]: (i) the nuclear incompressibility K with saturation density ρ c and the peak of the giant monopole resonance; (ii) the isoscalar effective mass M * s with the peak of the giant quadrupole resonance; (iii) the isovector effective mass M * v with the peak of the giant dipole resonance; and (iv) a block of correlated isovector indicators [39]: a sym , L, d ρn (E/N ), α D , and r n − r p . As seen in Fig. 5, unedf2 is missing some of the correlations predicted by SV-min. The reason is essentially that unedf2 has three symmetric neutron matter parameters fixed. Two of them, the isovector effective mass M * v and the slope of symmetry energy L, constitute crucial constraints because they are related to the properties of the linear response. Since they are fixed, they have been eliminated from the correlation matrix (white rows and columns). This step leaves the peak of the GDR unconstrained but constrains a sym considerably (uncertainty of 0.3 MeV for unedf2, compared with 1.7 MeV for SVmin). Consequently, nearly all correlations between a sym and the isovector static observables have disappeared. Another consequence of freezing L is the relatively small uncertainty for the neutron skin. These are due to the fact that the largest contribution to the error budget of r n comes from L [43].

Strength
The unedf2 value of a sym is consistent with the current estimates [46,76], and the same holds for α D and r n − r p [39,44,45,75]. A word of caution is in order concerning the peak position of the GDR. The energies of giant resonance peaks given in Table VII stem from an average over a broad energy region. Figure 6 shows the detailed energy-weighted dipole strength computed for 208 Pb with unedf2. The RPA results are folded with an energy-dependent width in order to simulate the increase of collisional width with excitation energy [44]. The GDR peak is strongly fragmented because it resides in a region of large density of 1ph states. The 1ph fragmentation is asymmetric because the density of 1ph states increases with energy. This produces a discrepancy between the averaged excitation energy E aver (the average taken just over the resonance region by virtue of a fluid dynamics approach) and the peak energy E peak , which is considerably smaller. The experimental energy of the GDR resonance, 13.6 MeV, corresponds to the peak energy. Consequently, we find that unedf2, similar to SV-min and several other Skyrme EDFs, underestimates the GDR peak energy. To overcome this problem, a smaller isovector effective mass or larger TRK sum rule enhancement is required [13,44].

C. Shell Structure
One of the primary motivations behind this work was to use experimental data on s.p. splittings to optimize the tensor coupling constants of the Skyrme EDF. In Table VIII, we report the root-mean-square deviations from experimental data for binding energies for 24 odd-A nuclei that are one mass unit away from the doubly magic systems 16 O, 40,48 Ca, 56 Ni, 132 Sn, and 208 Pb. (In the following, the abbreviation RMSD will always stand for a root-mean-square deviation between theoretical values and experimental data or empirical estimates.) The table also shows the RMSD for (six) twoneutron and two-proton separation energies across each shell gap. Since s.p. splittings are computed from binding energy differences of the neighboring odd-A nuclei, all these RMSDs are indicators of the quality of the underlying single-particle spectra. One can see in the table that although unedf0 and unedf1 reproduce binding energies B(A mag ± 1) better than SLy4, the latter works better for two-particle sepa-ration energies. The unedf2 parameterization brings a significant improvement on binding energies with respect to unedf0 while maintaining a decent reproduction of two-particle separation energies. In spite of this progress, the resulting RMSDs are quite appreciable.     48 Ca is clearly too small with unedf2. Otherwise, the positions of most of the levels seem to be slightly improved compared with unedf1, which was itself a minor improvement over un-edf0. The single-particle proton levels in 208 Pb show that the Z = 82 magic gap is also too small in unedf1 and unedf2, because of a low energy of the h 9/2 shell. Further, we notice in Fig. 9 the inversion of the 1j 15/2 and 1i 11/2 shells and a large shift in the energy of 3p 3/2 shell. The spectra shown in Figs. 7-9 are quite representative of the predictive power of the unedf family with respect to shell structure. To quantify further the quality of the predicted shell structure, we list in Table IX the RMSDs of singleparticle energies from the empirical values of Ref. [53]. The calculation is based on 75 (negative-energy) levels in the same set of double-magic nuclei as in Table VIII. We have also partitioned the set of nuclei into light (A < 80; 36 levels) and heavy nuclei (A ≥ 80; 39 levels). Note that all s.p. states used to compute the RMSDs were obtained from HFB calculations with the blocking procedure.
Overall, the RMSD from experiment is similar for all unedf parameterizations. The larger RMSD obtained for light nuclei is explained mostly by a lower level density, which increases the average error. Also, the impact of correlations missing in the Skyrme EFT approach is greater in lighter systems, the structure of which is profoundly impacted by surface effects. Even though twoparticle separation energies across the shell gap are improved with unedf2, the overall reproduction of shell structure is not.
These results are consistent with the conclusions of Ref. [26], where it was found that Skyrme EDFs are intrinsically limited in their ability to reproduce s.p. spectra in doubly-magic nuclei. The regression analysis technique employed therein suggests that the best possible RMSD for s.p. energies obtained in the Skyrme EDF approach is around 1.2 MeV. Although the calculations of Ref. [26] were performed at the HF level, it is unlikely that using the physically better motivated blocking procedure, and considering particle-vibration-coupling and self-interaction corrections [24] would significantly alter the conclusions. The RMSD of 1.38 MeV found for un-edf2 is thus very close to the limit given by the regression analysis, especially considering the diversity of constraints imposed during the fit.

D. Global Mass Table
The ability to reproduce nuclear properties globally across the whole nuclear landscape is one of the key requirements for an universal nuclear EDF. We have calculated the unedf2 nuclear mass table using the deformed HFB framework outlined in Ref. [77]. Figure 10 shows the residuals of the nuclear binding energies calculated with unedf2 with respect to the experimental values for isotopic and isotonic chains of even-even nuclei. Whereas the residuals for the isotopic chains show the typical arclike features common to many EDF calculations, these are hardly present in the isotonic chain residuals. It is difficult to explain this result, which may point to beyond mean-field effects not included in our functional and the related bias of the optimization [78]. Figure 11 shows the residuals obtained in unedf2 for two-neutron and two-proton separation energies. When compared with the prediction of unedf1 [15], the slightly worse RMSD reported in Table X primarily comes from larger deviations at the ends of each isotopic chain. As far as S 2p values are concerned, unedf1 yields values that are systematically too high. This trend is much less pronounced with unedf2. Table X lists the RMSDs for binding energies, twoparticle separation energies, pairing gaps, and proton radii of even-even nuclei. Compared with unedf1, un-edf2 is slightly less predictive for binding energies, S 2n values, and proton radii, but offers better reproduction of two-proton separation energies and neutron pairing gaps. The differences are, however, small.

E. Fission Barriers and Deformation Properties
One of the major differences between the original version of the unedf optimization protocol, used to determine the unedf0 parameterization, and its successive incarnations used to produce unedf1 and unedf2, is the inclusion of data on fission isomer excitation energies. This was motivated by the realization that surface properties of the energy density play a critical role in the EDF's ability to predict fission properties such as barriers and, consequently, spontaneous fission half-lives [79][80][81]. It was later shown that adding data corresponding to large nuclear deformations provides an effective constraint on the surface terms [82].
In Fig. 12, we present the residuals for the inner fission barrier heights, fission isomer excitation energies, and outer fission barrier heights in the actinide region calculated with unedf1, unedf2, the Gogny D1S model [81], and the Finite-Range Liquid Droplet Model (FRLDM) [85]. Although excitation energies of fission isomers are observables, fission barriers are not. Furthermore, the uncertainty on the empirical barrier heights ranges from 0.3 MeV [83] to 1 MeV, while the uncertainty for fission isomer energies ranges from 0.5 keV for 238 U to 0.5 MeV for 240 Pu (due to two different values reported in the literature) [84]. To keep the figure legible while conveying information on experimental uncertainties, the shaded area shows the average empirical error over the isotopes considered. All calculations were performed with the DFT solver hfodd of Ref. [86]. Details of the numerical implementation are discussed in Refs. [87,88]. As seen in Fig. 12, the deformation properties of the unedf2 functional are slightly degraded as compared to unedf1, especially for the outer barrier. The overall trend is that both barrier heights tend to be overestimated. This is quantified in Table XI, which lists the calculated RMSDs for the calculated first and second barrier heights, and fission isomer bandheads. The deviation from empirical values has increased by nearly 50% for the first barrier, and has doubled for the second barrier. The overall quality of unedf2 is now comparable to the SkM* parameterization [79].
As discussed in Ref. [82], the surface and surfacesymmetry coefficients of the leptodermous expansion of the nuclear energy determine average deformation properties of EDFs at large neutron-proton asymmetries. Table XII lists the coefficients of the liquid drop expansion extracted for the three unedf functionals, determined according to the methodology of Ref. [89]. We remark that the surface and curvature coefficients of both un-edf1 and unedf2 are very similar. However, the surfacesymmetry coefficient is significantly larger for the un-edf2 parametrization, and takes a value that is comparable to that of unedf0 and SkM*. This result explains why fission barriers (especially the outer barrier) are overestimated and similar to what can be obtained with SkM*.  fects and bulk properties that the EDF optimization has difficulties in keeping under control. As is well known, the spherical shell structure plays a major role in driving deformation properties [20]. Looking back at Fig. 9, we see that the positions of the neutron 1j 15/2 and proton 1i 13/2 shells in unedf2 are depleted as compared to experiment and unedf1. These high-j orbitals are especially sensitive to the surface terms of the functional and play an essential role in determining deformation properties of actinides.

F. Neutron Droplets
Trapped neutron droplets constitute a useful theoretical laboratory to test various many-body methods and effective interactions in inhomogeneous neutron matter. In particular, they probe the isovector channels of interactions or functionals, the role of which increases with neutron excess. The physics of neutron-rich nuclei is particularly relevant in the context of the inner crust of neutron stars [90], the r-process of nucleo-synthesis [90], and the determination of the limits of nuclear stability [16,41].
Since pure neutron matter is not self-bound, the neu-tron droplet must be confined by an external potential in order to produce bound states [91]. Recently, trapped neutron droplets have been used to test various ab initio approaches against DFT calculations with phenomenological functionals [92]. In particular, in Ref. [12], neutron droplets were used to test density matrix expansion techniques, which aim at building EDFs from the realistic interactions used in ab initio methods.
In Fig. 13, the binding energy per neutron of neutron droplets calculated with unedf0, unedf1, and un-edf2 are compared with the ab initio results obtained in Ref. [92] within the Auxiliary Field Diffusion Monte-Carlo (AFDMC) method. AFDMC calculations were performed with the AV8' parameterization of the twobody potential and the Urbana IX three-body interaction [93]. The figure also shows DFT calculations with SLy4, as well as a modified SLy4 parameterization that has been slightly readjusted in the isovector channel to reproduce the AFDMC results. All neutron droplet systems considered in Fig. 13 were confined by a spherical HO potential, with two choices of the oscillator frequency, ω = 5 MeV and ω = 10 MeV. As previously seen for unedf0 and unedf1 [15], unedf2 results are close to the ab initio calculations, even though the optimization did not include any information about neutron droplets. However, we notice that the results for N > 38 with ω = 10 MeV are not available for unedf2. This situation is the direct consequence of the neutron matter instabilities discussed in Sec. IV A. For N > 38 droplets, the central neutron density exceeds the critical density shown in Fig. 4; as a result, the HFB calculation fails to converge. For ω = 5 MeV, the central neutron density is low enough for higher particle numbers, so that the instabilities do not appear.

V. CONCLUSIONS
In this study, we have introduced the unedf2 parameterization of the Skyrme energy density. Compared with our previous work, there are two main differences: (i) we released the requirement that the isoscalar and isovector tensor coupling constants be zero, and (ii) we included experimental data on s.p. level splittings in doubly magic nuclei to better constrain spin-orbit and tensor coupling constants. In addition to these major changes, we have slightly extended our dataset to improve the pairing properties of the functional, especially in heavy nuclei. Following previous unedf optimizations, we have performed a comprehensive sensitivity analysis of our parameterization in order to obtain standard deviations and correlations among EDF parameters.
Global nuclear properties computed with unedf2 reflect little or no improvement with respect to our previous parameterizations. While the linear response analysis has shown that unedf2 does not have any finite-size instabilities in symmetric nuclear matter for densities up to 1.5ρ c , some instabilities are encountered in pure neutron matter, with the consequence that neutron droplet calculations do not converge at large neutron numbers and large oscillator frequencies. The position of the GDR peak in 208 Pb is slightly too low in energy, which is attributed to a persistent lack of constraints on the isovector effective mass. The quality of the single-particle shell structure near closed shell nuclei is almost as good as one can get with Skyrme EDFs, but this was almost the case with unedf0 and unedf1. The RMSD for nuclear binding energies is 1.95 MeV, which is far from the performance of semi-phenomenological mass models (see, for example, Ref. [94] for the most recent numbers) and comparable to unedf1. Deformation properties, which had been significantly improved with unedf1 are degraded markedly for unedf2, which yields fission barriers similar to that of the traditional SkM* functional.
On the other hand, as discussed in Sec. III D, the interval of confidence for the parameters is narrower for unedf2 than it was for unedf1, which itself was more tightly constrained than unedf0. In addition, the results of the sensitivity analysis of Sec. III E show that there is relatively weak dependence on individual experimental points. These results point to the fact that the coupling constants of the unedf2 functional are properly constrained by the data.
Although one can certainly improve the optimization protocol, for example by changing the relative weights in the χ 2 objective function, we believe this relative lack of improvement should be viewed as an intrinsic limitation of the Skyrme energy density, a local energy density that is up to second order in derivatives [48,49]. Indeed, as shown in Figs. 10-12, the residuals of various quantities predicted with unedf2 do not have a statistical distribution; hence, adding more data points or playing with the χ 2 is not going to change the situation as the deviations are mainly affected by systematic errors, i.e., imperfect modeling. In this context, unedf2 is an allaround Skyrme EDF that is fairly well constrained by various data, but it also marks the end of the Skyrme EDF strategy.
At this phase of nuclear DFT developments, it thus seems urgent to go beyond traditional Skyrme functionals. Two major avenues are being explored: one following the spirit of DFT, where the primary building block is the energy density functional that includes all correlation effects, and the other following the spirit of the selfconsistent mean-field theory, where the major ingredient is an effective pseudopotential and the beyond-mean-field correlations are added afterwards. The DFT description is especially convenient for tying in the energy density to a more fundamental theory of nuclear forces based, for example, on the chiral effective field theory. This can be accomplished by using EDF built from the density matrix expansion of realistic interactions [9][10][11][12]. A complementary route is to explore functionals with higher order derivatives of the density [95][96][97]. These EDFs are much richer than the Skyrme or Gogny functionals; hence, they should be able to capture more physics and reduce systematic errors.

ACKNOWLEDGMENTS
We are deeply indebted to the late M. Stoitsov, whose contribution to this work, especially the DFT solver and its interface with the pounders algorithm, was considerable. This work was supported by the U.S. Department of Energy under Contract Nos. de-sc0008499, DE-FG02-96ER40963, and DE-FG52-09NA29461 (University of Tennessee), No. DE-AC02-06CH11357 (Argonne National Laboratory), and DE-AC52-07NA27344 (Lawrence Livermore National Laboratory); by the Academy of Finland under the Centre of Excellence Programme 2012-2017 (Nuclear and Accelerator Based Physics Programme at JYFL) and FIDIPRO programme; and by the European Union's Seventh Framework Programme ENSAR (THEXO) under Grant No. 262010. Computational resources were provided through an INCITE award "Computational Nuclear Structure" by the National Center for Computational Sciences (NCCS) and National Institute for Computational Sciences (NICS) at Oak Ridge National Laboratory, through an award by the Livermore Computing Resource Center at Lawrence Livermore National Laboratory, and through an award by the Laboratory Computing Resource Center at Argonne National Laboratory.