Nuclear energy density optimization: Shell structure
Abstract
 Background

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.
 Purpose

In this work, we propose a new parameterization unedf2 of the Skyrme energy density functional.
 Methods

The functional optimization is carried out using the pounders optimization algorithm within the framework of the Skyrme HartreeFockBogoliubov 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 onebody density matrix. In order to impose constraints on all the parameters of the functional, selected data on singleparticle splittings in spherical doublymagic nuclei have been included into the experimental dataset.
 Results

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 singleparticle 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 spinorbit and tensor terms.
 Conclusions

unedf2 can be viewed as an allaround 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.
pacs:
21.10.k, 21.30.Fe, 21.60.Jz, 21.65.MnI Introduction
An important goal in research in lowenergy 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 SciDAC2 collaboration Bertsch et al. (2007); *(Fur11); *(Nam12); Bogner et al. (2013); its current successor, the NUCLEI SciDAC3 collaboration NUC (); and the FIDIPRO collaboration FID ().
Most of the nuclear EDFs used in selfconsistent meanfield calculations have been derived from phenomenological effective interactions, or pseudopotentials Vautherin and Brink (1972); Negele and Vautherin (1972); Bender et al. (2003). 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 Lalazissis et al. (2004); Gebremariam et al. (2010); Carlsson and Dobaczewski (2010); Stoitsov et al. (2010); Bogner et al. (2011). In a parallel effort, methodologies have been developed to validate nuclear EDFs by optimizing their lowenergy coupling constants to experimental data on finite nuclei and pseudodata on nuclear matter and other relevant systems Bogner et al. (2013); Klüpfel et al. (2009); Kortelainen et al. (2010a, 2012); Erler et al. (2013). 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 lowenergy 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 Bender et al. (2003); Stone and Reinhard (2007). 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 Bohr and Mottelson (1975). In an independentparticle picture, shell structure can be associated with the singleparticle (s.p.) spectra of the meanfield potential Ring and Schuck (2000); Nilsson and Ragnarsson (1995). 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 Pandharipande et al. (1997); Ring and Schuck (2000); Duguet and Hagen (2012) since s.p. motion is significantly modified by correlations Coló et al. (2010); Duguet and Hagen (2012); Tarpanov et al. (2014). 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 doublymagic nuclei Bender et al. (2003).
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, spinorbit, and tensor terms Satuła et al. (2008); Kortelainen et al. (2008). Suggestions to study tensor interactions within the selfconsistent meanfield approach were made already in the seventies Stancu et al. (1977) 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 Brown et al. (2006); Dobaczewski (2007); *(Dob07b); Lesinski et al. (2007); Zou et al. (2008); Tarpanov et al. (2008); Kortelainen et al. (2008); Satuła et al. (2008); Zalewski et al. (2008, 2009); Bender et al. (2009); MorenoTorres et al. (2010); Wang et al. (2011); Grasso and Anguiano (2013). 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 HartreeFockBogoliubov (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 Kortelainen et al. (2010a, 2012), 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. Kortelainen et al. (2010a), 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 oddeven 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. Kortelainen et al. (2012), we modified the form of the functional by removing the center of mass correction, which allows for straightforward timedependent HartreeFock (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 onebody 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 doublymagic 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 Reinhard and Nazarewicz (2010); Fattoyev and Piekarewicz (2011); Erler et al. (2012a); Gao et al. (2013); Kortelainen et al. (2013); Reinhard and Nazarewicz (2013); Reinhard et al. (2013); Nazarewicz et al. (2014). Since the previous work of Ref. Kortelainen et al. (2008) based on the linear regression methodology demonstrated that the current standard form of the Skyrme EDF cannot ensure a spectroscopicquality description of s.p. energies, unedf2 is certainly not the universal nuclear EDF. It can be viewed, however, as the best allaround 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 is a functional of the onebody density matrix and pairing density and can be cast into the generic form
(1)  
where is the kinetic energy; is the isoscalar () and isovector () particlehole Skyrme energy density; is the pairing energy density; and is the Coulomb term.
The particlehole part of the Skyrme energy density reads
(2)  
where each term is multiplied by a coupling constant represented by a real number. The coupling constant is the only exception, as it has the traditional densitydependence
(3) 
The definitions of the various densities , , and ( is the vector part of ) can be found in Ref. Engel et al. (1975); Dobaczewski and Dudek (1996); Bender et al. (2003); Perlińska et al. (2004); Lesinski et al. (2007). The coupling constants and are related to the volume part of the energy density and can be expressed as a function of the parameters of infinite nuclear matter Kortelainen et al. (2010a).
The term in Eq. (2) represents a tensor energy density. The approximation that was made in the unedf0 and unedf1 optimizations was to set for both and . 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 : each of these components is multiplied by the same coupling constant, see Dobaczewski and Dudek (1996); Perlińska et al. (2004).
The pairing term is derived from the mixed pairing force of Ref. Dobaczewski et al. (2002), leading to the energy density
(4) 
where is the pairing strength. In this work as before, we take =0.16 fm. We have used different pairing strengths for neutrons and protons Bertsch et al. (2009). Owing to the zero range of our effective pairing force, we have used a pairing cutoff MeV to truncate the quasiparticle space. To prevent the collapse of pairing correlations near closed shells, we have also used the variant of the LipkinNogami (LN) method as in Ref. Stoitsov et al. (2003).
As with unedf1, we disregard the centerofmass correction; for motivation see Ref. Kortelainen et al. (2012). The Coulomb term contains direct and exchange contributions. The direct part is computed from the proton density distribution assuming a pointproton 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.1, we review the experimental data used in the fit. In Sec. III.2, we describe how information on shell structure was incorporated. The optimization procedure itself is summarized in Sec. III.3. Section III.4 compares unedf2 to previous unedf parameterizations, and Sec. III.5 contains the results of the sensitivity analysis.
iii.1 Experimental Dataset
In order to determine the tensor coupling constants , 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 singleparticle level splittings, five new data points for oddeven staggering (OES), and one additional binding energy.
Nucleus  Level  Energy  

Ca  6.80  
Ca  7.28  
Ca  7.24  
Ca  8.80  
Ca  4.92  
Sn  6.68  
Sn  6.03  
Pb  6.08  
Pb  5.56 
We show in Table 1 the empirical values of singleparticle splittings in several doublymagic nuclei. All values are taken from the empirical s.p. energies listed in Ref. Schwierz et al. (2007) except for the proton splitting in Ca, which is taken from Ref. Oros (1996). 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 singleparticle data points in the function to MeV. This choice was motivated based on the singular value decomposition (SVD) analysis performed in Ref. Kortelainen et al. (2008), 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 Sn requires the groundstate energy of Sn, see Sec. III.2. For this reason, we added the binding energy of this nucleus to the dataset. As in Ref. Kortelainen et al. (2010a), 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 MeV. For this additional datum, we took the same weight MeV as for other binding energies.
Neutrons  Protons  

90  142  0.681450  90  142  0.813287 
50  74  1.250400  76  90  1.169046 
50  70  1.316825 
In addition to the singleparticle splittings and the binding energy of Sn, we have added five new OES data points, which are listed in Table 2. This was motivated by the observation that pairing properties of actinide nuclei and neutronrich tin isotopes are poorly reproduced by unedf1, suggesting that the weight of pairingrelated data in the objective function should be increased. We recall that the experimental OES that we use is defined as the average of two oddeven (protons) and evenodd (neutrons) values, that is, . In addition to these new experimental points, we have increased the weight of all OES data points in the optimization from to .
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 singleparticle level splittings. The changes with respect to the unedf1 optimization procedure are as follows:

Tensor coupling constants and are optimized;

The binding energy of Sn is added to dataset;

9 new singleparticle splittings are included in dataset;

5 new OES data points are added to dataset;

The weight of all OES data points is increased to .
iii.2 Computation of SingleParticle 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. Alex Brown (1998); Chabanat et al. (1998); Klüpfel et al. (2009)). Within this approximation, the manybody wavefunction reduces to a single Slater determinant, and the theoretical singleparticle levels are taken as the eigenvalues of the HF Hamiltonian following Koopmans’ theorem Koopmans (1934); Kortelainen et al. (2008). In the framework of the HFB theory with approximate particle number projection, where the basic degrees of freedom are not particles but quasiparticles, separation energies provide a more convenient quantity to relate to effective singleparticle 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 equalfilling approximation to blocking, since it yields an excellent estimate of the full symmetrybreaking blocking results Schunck et al. (2010). Since the shape polarization induced by the blocking prescription spontaneously breaks the spherical symmetry of the odd nucleus, the energy degeneracy of a blocked spherical quasiparticle orbital with angular momentum is lifted Zalewski et al. (2008); Schunck et al. (2010); Satuła et al. (2012). In the equal filling approximation, however, timereversal symmetry is conserved, and states with the angular momentum projection + and  are degenerate.
This fragmentation of any given 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 can originate from a 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 . The associated systematic error does not exceed 100 keV Kortelainen et al. (2012).
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 oneparticle 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 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 state in Ca is due to the strong fragmentation of the strength among multiple states, resulting in a rather broad centroid Uozumi et al. (1994). Similarly, the proton state in Ca is fragmented among multiple states Burrows (2008). In the Supplemental Material at SM (), we provide both sets of experimental s.p. splittings for convenience.
Our procedure to generate theoretical s.p. splittings follows that of Ref. Rutz et al. (1998) 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 and perform the blocking calculation in the system with particles. The effective particle and hole s.p. energies are respectively defined as
(5a)  
(5b) 
where refers to the particle number of the reference (doublymagic) nucleus and 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 doublymagic 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 eveneven binding energy cancels out, and the s.p. splitting reduces to the difference of total binding energies of blocked configurations.
iii.3 Optimization
All HFB calculations in the unedf2 optimization were performed with the DFT solver hfbtho Stoitsov et al. (2013). 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 superdeformed (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 . In this work, we have maintained the same setup as for unedf1: the spherical basis is used for all groundstate configurations, and the stretched basis with is used for SD states. In all cases, the spherical frequency of the HO basis is set at Stoitsov et al. (2013).
The objective function in our optimization is
(6) 
where is the number of different data types; for unedf2 . The total number of data points is (here ), and the number of parameters to be fitted is (here ). The calculated value of the observable of type is , while the corresponding experimental value is denoted by . Each data type has a weight . As in the cases of unedf0 and unedf1, the isovector effective mass 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 for historical reasons Chabanat et al. (1998). In the Supplemental Material SM (), 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 . During the optimization, two of the parameters, and , 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. Kortelainen et al. (2010a) for details.
iii.4 The UNEDF2 parameterization
The optimized parameter set of the EDF unedf2 is listed in Table 3 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 SM () 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.
95% CI  

0.15631  0.00112  [ 0.154, 0.158]  
15.8      
239.930  10.119  [ 223.196, 256.663]  
29.131  0.321  [ 28.600, 29.662]  
40.0      
1.074  0.052  [ 0.988, 1.159]  
46.831  2.689  [ 51.277, 42.385]  
113.164  24.322  [153.383, 72.944]  
208.889  8.353  [222.701,195.077]  
230.330  6.792  [241.561,219.099]  
64.309  5.841  [ 73.968, 54.649]  
38.650  15.479  [ 64.246, 13.054]  
54.433  16.481  [ 81.687, 27.180]  
65.903  17.798  [ 95.334, 36.472] 
We first note that for unedf2, the nuclear incompressibility parameter , while in the top range of acceptable values, is now constrained by the data, whereas the slope of the symmetry energy 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 4 lists all three unedf parameterizations produced so far, and compares them to the SLy4 parametrization, which was the starting point for unedf0. 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 spinorbit coupling constants. Indeed, it was shown in Ref. Lesinski et al. (2007) that there is a strong anticorrelation between the isoscalar spinorbit and tensor coupling constants. This relationship is confirmed in our optimization through a large correlation coefficient of between and . In fact, the values of and are consistent with the empirical dependence reported in Ref. Lesinski et al. (2007). Yet, in spite of this very strong correlation, the value of changes only by 13% between unedf1 and unedf2.
SLy4  unedf0  unedf1  unedf2  

0.16000  0.16053  0.15871  0.15631  
15.972  16.056  15.8  15.8  
229.901  230.0  220.0  239.930  
32.004  30.543  28.987  29.131  
45.962  45.080  40.005  40.0  
1.439  0.9  0.992  1.074  
76.996  55.261  45.135  46.831  
+15.657  55.623  145.382  113.164  
258.200  170.374  186.065  208.889  
258.200  199.202  206.580  230.330  
92.250  79.531  74.026  64.309  
30.750  45.630  35.658  38.650  
0.000  0.000  0.000  54.433  
0.000  0.000  0.000  65.903 
Looking more closely at the values of all spinorbit and tensor coupling constants of the unedf2 parameterization, we find that they are compatible with the results of Ref. Zalewski et al. (2008). In particular, the value of is close to the “universal” value of MeV obtained from a refit of three different EDFs to the spinorbit splittings in Ca and Ni. Our optimized tensor coupling constants, as well as their isoscalarisovector trend, are also in the same ballpark as those partial refits of Ref. Zalewski et al. (2008). They are, however, just beyond the space of the T family of parameterizations considered in Ref. Lesinski et al. (2007). Recent work also suggests that only the region where and should be physical Grasso and Anguiano (2013). 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 spinorbit and tensor terms.
Channel  

0.733  0.791  0.134  0.639  0.878  0.743  
0.328  0.291  0.319  1.545  0.528  0.900 
We show in Table 5 the unedf2 coupling constants in natural units Furnstahl and Hackworth (1997); Kortelainen et al. (2010b). 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 , which was found in Ref. Kortelainen et al. (2010b) to be valid for Skyrme EDFs. As seen in Table 5, nearly all the unedf2 coupling constants fall in the interval which is compatible with the hypothesis of naturalness Kortelainen et al. (2010b). The one notable exception is , which is unnaturally small; and are also at the limits of the allowed interval.
iii.5 Sensitivity Analysis
The standard deviations of the unedf2 parameterization are listed in Table 3, 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.
1.00  
0.97  1.00  
0.07  0.03  1.00  
0.08  0.05  0.24  1.00  
0.43  0.43  0.22  0.89  1.00  
0.42  0.37  0.83  0.17  0.31  1.00  
0.06  0.02  0.27  0.96  0.85  0.17  1.00  
0.09  0.05  0.21  0.89  0.80  0.14  0.86  1.00  
0.51  0.50  0.34  0.40  0.68  0.55  0.36  0.34  1.00  
0.31  0.29  0.19  0.00  0.04  0.18  0.07  0.02  0.14  1.00  
0.56  0.55  0.26  0.05  0.35  0.53  0.02  0.02  0.88  0.35  1.00  
0.36  0.35  0.13  0.23  0.16  0.14  0.29  0.25  0.02  0.57  0.29  1.00  
Table 6 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 , 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 and . This is reminiscent of our first parameterization, dubbed unedfnb in Ref. Kortelainen et al. (2010a), where no bounds had been imposed on the coupling constants. Further, as discussed in Sec. III.4, the isoscalar spinorbit and tensor coupling constants are strongly anticorrelated and seem to follow the trend predicted in Ref. Lesinski et al. (2007).
The overall impact of each data type on the unedf2 parameterization can be assessed by studying the sensitivity matrix :
(7) 
where is the Jacobian matrix calculated with the parameterization . Following Refs. Kortelainen et al. (2010a, 2012), we have calculated the partial sums of the absolute values in 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 . 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., , , and ) 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 (, , ) behave differently but are less constrained by the data, as shown by their large standard deviations reported in Table 3.
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
(8) 
for the optimal solution is presented when each data point is shifted by an amount of . The standard deviations for parameters are listed in Table 3. As for unedf0 and unedf1, the overall changes in are of the same order of magnitude and are very small, . 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 Ca.
Iv Properties of UNEDF2 parameterization
In this section, we review various properties of the unedf2 parameterization. In Sec. IV.1, we apply the linear response theory to test the functional against the presence of finitesize instabilities. Section IV.2 discusses correlations among various observables. Predictions of unedf2 for shell structure in doublymagic nuclei are presented in Sec. IV.3, and global binding energy and deformation trends (in particular in the context of nuclear fission) are analyzed in Secs. IV.4 and IV.5, respectively. Section IV.6 contains the discussion of unedf2 predictions for neutron droplets in external traps.
iv.1 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. Pastore et al. (2012a); *(Pas12b) 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. Perlińska et al. (2004). This recent development is needed here because of the nonstandard spinorbit term of the unedf parametrizations. The response functions of interest to the present work are defined as the response of the infinite medium to external probes of the type
(9) 
where () is the spin (its projection along the axis), is the isospin, and is an operator acting on spin, isospin, or both. As recently shown in Ref. Pastore et al. (2012c), the response function of neutron matter can provide information about instabilities in finite nuclei Lesinski et al. (2006). More precisely, it was shown that whenever a pole appears in the response function close enough to the saturation density of the system , the finite nucleus undergoes an instability in the corresponding channel (scalar/vector, isoscalar/isovector), see examples in Ref. Schunck et al. (2010).
Several quantitative criteria to estimate the likelihood of finitesize instabilities for a given EDF have been recently proposed Pastore et al. (2013); Hellemans et al. (2013). Because of shell fluctuations, the nucleus can explore regions of densities slightly larger than the saturation density. In Ref. Hellemans et al. (2013), the following conservative criterion was established: whenever the response function has a pole at a density , there is a risk of instability in calculations for finite nuclei. More complex criteria were proposed in Pastore et al. (2012b). The poles of the response function are in practice determined by solving the equation
(10) 
where the expressions of can be found in Ref. Pastore et al. (2012a).
Figure 3 shows the position of the lowest poles of the response function as a function of the transferred momentum in symmetric nuclear matter for unedf0, unedf1, unedf2, and SkP Dobaczewski et al. (1984). Since the unedf functionals have been developed to be used only in the timeeven channel, we will limit our analysis to . We first remark that in all cases we have an instability in the region at low density and low momentum in the channel . This is the wellknown spinodal instability, which is physical. We then observe that all unedf functionals satisfy the stability criterion given in Ref. Hellemans et al. (2013). On the other hand, SkP is unstable in the scalar/isovector channel in the region of densities around ; as a result it was shown in Ref. Lesinski et al. (2006) 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 instabilities could manifest in neutronrich nuclei or in trapped neutron droplets. As we will see in Sec. IV.6 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 Lesinski et al. (2007).
iv.2 Correlations with Other Observables
The sensitivity analysis presented in Sec. III.5 aimed at quantifying the behavior of the 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 Reinhard and Nazarewicz (2010); Gao et al. (2013); Kortelainen et al. (2013). 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 Brandt (1999). It is a slightly different prescription from the procedure that we use to define the confidence intervals; see Sec. III.B.1. of Ref. Kortelainen et al. (2010a) for details. Asymptotically (at the limit of ), both prescriptions are in fact equivalent (see discussion in Sec. 3.3.1 of Ref. Seber and Wild (1989)).
Observable  

0.156  0.001  
240  10  
0.93  0.04  
29.1  0.3  
75.8  2.2  
GMR  14.0  0.3 
GQR  10.8  0.3 
GDR  13.6  0.1 
13.8  0.1  
0.167  0.003 
Here we construct an approximate confidence region using the following approach. At the minimum , the quantity characterizes the bestfit parameterization. The parameters in a neighborhood of the minimum can still provide a reasonable description of nuclear properties. The parameter space is thus referred to as a “reasonable” domain. Each observable that can be computed in the EDF theory is also a function of the Skyrme parameters, . Varying in the vicinity of the optimal set will lead to fluctuations in the values of with respect to its value at the minimum, . The uncertainty of the prediction is characterized by the variance , where the average value is computed from
(11) 
This simple estimate of uncertainties provides valuable information on the predictive power of the model. Further information can be obtained from the correlation coefficient between two observables and defined from the covariance matrix as
(12) 
Table 7 shows the values and the uncertainties of a large set of observables. The vicinity was defined by the level set . As noted earlier, this construction of is different from the one used to define the 95% confidence interval, and hence, the standard deviations reported in Table 7 are slightly different from those of Table 3. As discussed later in Sec. III.5, the set unedf2 delivers rather small uncertainties for all observables shown.
Figure 5 shows the correlation matrix (12) between various pairs of observables computed for Pb. To illustrate the impact of the optimization protocol, we compare unedf2 with the SVmin parameterization Klüpfel et al. (2009). In the case of SVmin (upper panel), we can see four blocks of highly correlated observables Reinhard and Nazarewicz (2013); Piekarewicz et al. (2012); Nazarewicz et al. (2014): (i) the nuclear incompressibility with saturation density and the peak of the giant monopole resonance; (ii) the isoscalar effective mass with the peak of the giant quadrupole resonance; (iii) the isovector effective mass with the peak of the giant dipole resonance; and (iv) a block of correlated isovector indicators Reinhard and Nazarewicz (2010): , , , , and .
As seen in Fig. 5, unedf2 is missing some of the correlations predicted by SVmin. The reason is essentially that unedf2 has three symmetric neutron matter parameters fixed. Two of them, the isovector effective mass and the slope of symmetry energy , 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 considerably (uncertainty of 0.3 MeV for unedf2, compared with 1.7 MeV for SVmin). Consequently, nearly all correlations between and the isovector static observables have disappeared. Another consequence of freezing is the relatively small uncertainty for the neutron skin. These are due to the fact that the largest contribution to the error budget of comes from Kortelainen et al. (2013).
The unedf2 value of is consistent with the current estimates Lattimer (2012); Nazarewicz et al. (2014), and the same holds for and Reinhard and Nazarewicz (2010); Piekarewicz et al. (2012); Reinhard and Nazarewicz (2013); Reinhard et al. (2013). A word of caution is in order concerning the peak position of the GDR. The energies of giant resonance peaks given in Table 7 stem from an average over a broad energy region. Figure 6 shows the detailed energyweighted dipole strength computed for Pb with unedf2. The RPA results are folded with an energydependent width in order to simulate the increase of collisional width with excitation energy Reinhard and Nazarewicz (2013). The GDR peak is strongly fragmented because it resides in a region of large density of states. The fragmentation is asymmetric because the density of states increases with energy. This produces a discrepancy between the averaged excitation energy (the average taken just over the resonance region by virtue of a fluid dynamics approach) and the peak energy , 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 SVmin 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 Klüpfel et al. (2009); Reinhard and Nazarewicz (2013).
iv.3 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 8, we report the rootmeansquare deviations from experimental data for binding energies for 24 odd nuclei that are one mass unit away from the doubly magic systems O, Ca, Ni, Sn, and Pb. (In the following, the abbreviation RMSD will always stand for a rootmeansquare deviation between theoretical values and experimental data or empirical estimates.) The table also shows the RMSD for (six) twoneutron and twoproton separation energies across each shell gap. For example, in the case of Pb: would stand for (protons) and (neutrons). Similarly, for protons and for neutrons; represents the twoneutron separation energy of Pb: ; and is the twoproton separation energy of Bi . Since s.p. splittings are computed from binding energy differences of the neighboring odd nuclei, all these RMSDs are indicators of the quality of the underlying singleparticle spectra.
Quantity  SLy4  unedf0  unedf1  unedf2 

3.30  2.70  2.78  2.16  
3.06  2.44  2.04  2.12  
and  0.90  1.44  1.59  0.90 
One can see in the table that although unedf0 and unedf1 reproduce binding energies better than SLy4, the latter works better for twoparticle separation energies. The unedf2 parameterization brings a significant improvement on binding energies with respect to unedf0 while maintaining a decent reproduction of twoparticle separation energies. In spite of this progress, the resulting RMSDs are quite appreciable.
Figures 7–9 display the s.p. levels, as defined by Eqs.(5a)(5b), for neutrons in Ca and for protons and neutrons in Pb, respectively. Compared with the empirical values, the gap in 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 unedf0. The singleparticle proton levels in Pb show that the magic gap is also too small in unedf1 and unedf2, because of a low energy of the shell. Further, we notice in Fig. 9 the inversion of the 1 and 1 shells and a large shift in the energy of 3 shell. The spectra shown in Figs. 79 are quite representative of the predictive power of the unedf family with respect to shell structure.
Nuclei  unedf0  unedf1  unedf2 

All  1.42  1.38  1.38 
Light  1.80  1.72  1.74 
Heavy  0.94  0.97  0.95 
To quantify further the quality of the predicted shell structure, we list in Table 9 the RMSDs of singleparticle energies from the empirical values of Ref. Schwierz et al. (2007). The calculation is based on 75 (negativeenergy) levels in the same set of doublemagic nuclei as in Table 8. We have also partitioned the set of nuclei into light (; 36 levels) and heavy nuclei (; 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. Kortelainen et al. (2008), where it was found that Skyrme EDFs are intrinsically limited in their ability to reproduce s.p. spectra in doublymagic 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. Kortelainen et al. (2008) were performed at the HF level, it is unlikely that using the physically better motivated blocking procedure, and considering particlevibrationcoupling and selfinteraction corrections Tarpanov et al. (2014) would significantly alter the conclusions. The RMSD of 1.38 MeV found for unedf2 is thus very close to the limit given by the regression analysis, especially considering the diversity of constraints imposed during the fit.
iv.4 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. Erler et al. (2012b). 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 eveneven 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 meanfield effects not included in our functional and the related bias of the optimization Bender et al. (2006).
Figure 11 shows the residuals obtained in unedf2 for twoneutron and twoproton separation energies. When compared with the prediction of unedf1 Kortelainen et al. (2012), the slightly worse RMSD reported in Table 10 primarily comes from larger deviations at the ends of each isotopic chain. As far as values are concerned, unedf1 yields values that are systematically too high. This trend is much less pronounced with unedf2.
Table 10 lists the RMSDs for binding energies, twoparticle separation energies, pairing gaps, and proton radii of eveneven nuclei. Compared with unedf1, unedf2 is slightly less predictive for binding energies, values, and proton radii, but offers better reproduction of twoproton separation energies and neutron pairing gaps. The differences are, however, small.
Observable  unedf0  unedf1  unedf2  No. 

1.428  1.912  1.950  555  
2.092  2.566  2.475  113  
1.200  1.705  1.792  442  
0.758  0.752  0.843  500  
1.447  1.161  1.243  99  
0.446  0.609  0.711  401  
0.862  0.791  0.778  477  
1.496  1.264  1.309  96  
0.605  0.618  0.572  381  
0.355  0.358  0.285  442  
0.401  0.388  0.327  89  
0.342  0.350  0.273  353  
0.258  0.261  0.276  395  
0.346  0.304  0.472  83  
0.229  0.248  0.194  312  
0.017  0.017  0.018  49  
0.022  0.019  0.020  16  
0.013  0.015  0.017  33 
iv.5 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 halflives Bartel et al. (1982); Tondeur (1985); Berger et al. (1989). It was later shown that adding data corresponding to large nuclear deformations provides an effective constraint on the surface terms Nikolov et al. (2011).
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 Berger et al. (1989), and the FiniteRange Liquid Droplet Model (FRLDM) Möller et al. (2009). Although excitation energies of fission isomers are observables, fission barriers are not. Furthermore, the uncertainty on the empirical barrier heights ranges from MeV Smirenkin (1993) to MeV, while the uncertainty for fission isomer energies ranges from keV for U to MeV for Pu (due to two different values reported in the literature) Singh et al. (2002). 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. Schunck et al. (2012). Details of the numerical implementation are discussed in Refs. McDonnell et al. (2013); Schunck (2013).
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 11, 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 Bartel et al. (1982).
unedf2  unedf1  FRLDM  SkM*  D1S  

1.470  1.030  1.520  1.610  0.709  
0.515  0.357  0.675  0.351  0.339  
1.390  0.690  1.130  1.390  1.140 
As discussed in Ref. Nikolov et al. (2011), the surface and surfacesymmetry coefficients of the leptodermous expansion of the nuclear energy determine average deformation properties of EDFs at large neutronproton asymmetries. Table 12 lists the coefficients of the liquid drop expansion extracted for the three unedf functionals, determined according to the methodology of Ref. Reinhard et al. (2006). We remark that the surface and curvature coefficients of both unedf1 and unedf2 are very similar. However, the surfacesymmetry coefficient is significantly larger for the unedf2 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*.
Functional  

unedf0  16.056  30.543  18.7  7.1  44 
unedf1  15.800  29.987  16.7  8.8  29 
unedf2  15.800  29.131  16.8  8.7  42 
SkM*  15.752  30.040  17.6  9.0  52 
It also suggests a complex interplay between shell effects 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 Nilsson and Ragnarsson (1995). Looking back at Fig. 9, we see that the positions of the neutron 1 and proton 1 shells in unedf2 are depleted as compared to experiment and unedf1. These high orbitals are especially sensitive to the surface terms of the functional and play an essential role in determining deformation properties of actinides.
iv.6 Neutron Droplets
Trapped neutron droplets constitute a useful theoretical laboratory to test various manybody 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 neutronrich nuclei is particularly relevant in the context of the inner crust of neutron stars Ravenhall et al. (1983), the process of nucleosynthesis Ravenhall et al. (1983), and the determination of the limits of nuclear stability Erler et al. (2012a, 2013).
Since pure neutron matter is not selfbound, the neutron droplet must be confined by an external potential in order to produce bound states Pieper (2003). Recently, trapped neutron droplets have been used to test various ab initio approaches against DFT calculations with phenomenological functionals Gandolfi et al. (2011). In particular, in Ref. Bogner et al. (2011), 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 unedf2 are compared with the ab initio results obtained in Ref. Gandolfi et al. (2011) within the Auxiliary Field Diffusion MonteCarlo (AFDMC) method. AFDMC calculations were performed with the AV8’ parameterization of the twobody potential and the Urbana IX threebody interaction Pudliner et al. (1997); *(Pud96). 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, MeV and MeV. As previously seen for unedf0 and unedf1 Kortelainen et al. (2012), 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 with MeV are not available for unedf2. This situation is the direct consequence of the neutron matter instabilities discussed in Sec. IV.1. For droplets, the central neutron density exceeds the critical density shown in Fig. 4; as a result, the HFB calculation fails to converge. For 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 spinorbit 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 finitesize instabilities in symmetric nuclear matter for densities up to , 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 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 singleparticle 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 semiphenomenological mass models (see, for example, Ref. Goriely et al. (2013) 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.4, 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.5 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 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 Dobaczewski and Dudek (1996); Perlińska et al. (2004). Indeed, as shown in Figs. 1012, the residuals of various quantities predicted with unedf2 do not have a statistical distribution; hence, adding more data points or playing with the 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 meanfield theory, where the major ingredient is an effective pseudopotential and the beyondmeanfield 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 Stoitsov et al. (2010); Bogner et al. (2011); Gebremariam et al. (2010); Carlsson and Dobaczewski (2010). A complementary route is to explore functionals with higher order derivatives of the density Carlsson et al. (2008); Raimondi et al. (2011); Dobaczewski et al. (2012). These EDFs are much richer than the Skyrme or Gogny functionals; hence, they should be able to capture more physics and reduce systematic errors.
Acknowledgements.
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. DESC0008499, DEFG0296ER40963, and DEFG5209NA29461 (University of Tennessee), No. DEAC0206CH11357 (Argonne National Laboratory), and DEAC5207NA27344 (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.References
 Bertsch et al. (2007) G. Bertsch, D. Dean, and W. Nazarewicz, SciDAC Review 6, 42 (2007).
 Furnstahl (2011) R. Furnstahl, Nucl. Phys. News 21, 18 (2011).
 Nam et al. (2012) H. Nam, M. Stoitsov, W. Nazarewicz, A. Bulgac, G. Hagen, M. Kortelainen, P. Maris, J. C. Pei, K. J. Roche, N. Schunck, I. Thompson, J. P. Vary, and S. M. Wild, J. Phys.: Conf. Ser. 402, 012033 (2012).
 Bogner et al. (2013) S. Bogner, A. Bulgac, J. Carlson, J. Engel, G. Fann, R. Furnstahl, S. Gandolfi, G. Hagen, M. Horoi, C. Johnson, M. Kortelainen, E. Lusk, P. Maris, H. Nam, P. Navratil, W. Nazarewicz, E. Ng, G. Nobre, E. Ormand, T. Papenbrock, J. Pei, S. Pieper, S. Quaglioni, K. Roche, J. Sarich, N. Schunck, M. Sosonkina, J. Terasaki, I. Thompson, J. Vary, and S. Wild, Comput. Phys. Comm. 184, 2235 (2013).
 (5) http://computingnuclei.org.
 (6) https://www.jyu.fi/fysiikka/en/research/accelerator/fidipro.
 Vautherin and Brink (1972) D. Vautherin and D. M. Brink, Phys. Rev. C 5, 626 (1972).
 Negele and Vautherin (1972) J. W. Negele and D. Vautherin, Phys. Rev. C 5, 1472 (1972).
 Bender et al. (2003) M. Bender, P.H. Heenen, and P.G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
 Lalazissis et al. (2004) G. A. Lalazissis, P. Ring, and D. Vretenar, Extended Density Functionals in Nuclear Structure Physics (Lecture Notes in Physics 641, Springer, 2004).
 Gebremariam et al. (2010) B. Gebremariam, T. Duguet, and S. K. Bogner, Phys. Rev. C 82, 014305 (2010).
 Carlsson and Dobaczewski (2010) B. G. Carlsson and J. Dobaczewski, Phys. Rev. Lett. 105, 122501 (2010).
 Stoitsov et al. (2010) M. Stoitsov, M. Kortelainen, S. K. Bogner, T. Duguet, R. J. Furnstahl, B. Gebremariam, and N. Schunck, Phys. Rev. C 82, 054307 (2010).
 Bogner et al. (2011) S. K. Bogner, R. J. Furnstahl, H. Hergert, M. Kortelainen, P. Maris, M. Stoitsov, and J. P. Vary, Phys. Rev. C 84, 044306 (2011).
 Klüpfel et al. (2009) P. Klüpfel, P.G. Reinhard, T. J. Bürvenich, and J. A. Maruhn, Phys. Rev. C 79, 034310 (2009).
 Kortelainen et al. (2010a) M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild, Phys. Rev. C 82, 024313 (2010a).
 Kortelainen et al. (2012) M. Kortelainen, J. McDonnell, W. Nazarewicz, P.G. Reinhard, J. Sarich, N. Schunck, M. V. Stoitsov, and S. M. Wild, Phys. Rev. C 85, 024304 (2012).
 Erler et al. (2013) J. Erler, C. J. Horowitz, W. Nazarewicz, M. Rafalski, and P.G. Reinhard, Phys. Rev. C 87, 044320 (2013).
 Stone and Reinhard (2007) J. Stone and P.G. Reinhard, Prog. Part. Nucl. Phys. 58, 587 (2007).
 Bohr and Mottelson (1975) A. Bohr and B. Mottelson, Nuclear Structure, vol. II (W A. Benjamin, Reading, 1975).
 Ring and Schuck (2000) P. Ring and P. Schuck, The Nuclear ManyBody Problem (SpringerVerlag, 2000).
 Nilsson and Ragnarsson (1995) S. Nilsson and I. Ragnarsson, Shapes and Shells in Nuclear Structure (Cambridge University Press, Cambridge, 1995).
 Pandharipande et al. (1997) V. R. Pandharipande, I. Sick, and P. K. A. Huberts, Rev. Mod. Phys. 69, 981 (1997).
 Duguet and Hagen (2012) T. Duguet and G. Hagen, Phys. Rev. C 85, 034330 (2012).
 Coló et al. (2010) G. Coló, H. Sagawa, and P. Bortignon, Phys. Rev. C 82 (2010).
 Tarpanov et al. (2014) D. Tarpanov, J. Toivanen, J. Dobaczewski, and B. Carlsson, Phys. Rev. C 89, 014307 (2014).
 Satuła et al. (2008) W. Satuła, R. A. Wyss, and M. Zalewski, Phys. Rev. C 78, 011302 (2008).
 Kortelainen et al. (2008) M. Kortelainen, J. Dobaczewski, K. Mizuyama, and J. Toivanen, Phys. Rev. C 77, 064307 (2008).
 Stancu et al. (1977) F. Stancu, D. Brink, and H. Flocard, Phys. Lett. B 68, 108 (1977).
 Brown et al. (2006) B. A. Brown, T. Duguet, T. Otsuka, D. Abe, and T. Suzuki, Phys. Rev. C 74, 061303 (2006).
 Dobaczewski (2007) J. Dobaczewski, in Opportunities with Exotic Beams, edited by T. Duguet, H. Esbensen, K. Nollett, and C. Roberts (World Scientific, Singapore, 2007) p. 152.
 Dobaczewski et al. (2007) J. Dobaczewski, N. Michel, W. Nazarewicz, M. Płoszajczak, and J. Rotureau, Prog. Part. Nucl. Phys. 59, 432 (2007).
 Lesinski et al. (2007) T. Lesinski, M. Bender, K. Bennaceur, T. Duguet, and J. Meyer, Phys. Rev. C 76, 014312 (2007).
 Zou et al. (2008) W. Zou, G. Colò, Z. Ma, H. Sagawa, and P. F. Bortignon, Phys. Rev. C 77, 014314 (2008).
 Tarpanov et al. (2008) D. Tarpanov, H. Liang, N. V. Giai, and C. Stoyanov, Phys. Rev. C 77, 054316 (2008).
 Zalewski et al. (2008) M. Zalewski, J. Dobaczewski, W. Satuła, and T. R. Werner, Phys. Rev. C 77, 024316 (2008).
 Zalewski et al. (2009) M. Zalewski, P. Olbratowski, M. Rafalski, W. Satuła, T. R. Werner, and R. A. Wyss, Phys. Rev. C 80, 064307 (2009).
 Bender et al. (2009) M. Bender, K. Bennaceur, T. Duguet, P. H. Heenen, T. Lesinski, and J. Meyer, Phys. Rev. C 80, 064302 (2009).
 MorenoTorres et al. (2010) M. MorenoTorres, M. Grasso, H. Liang, V. De Donno, M. Anguiano, and N. Van Giai, Phys. Rev. C 81, 064327 (2010).
 Wang et al. (2011) Y. Z. Wang, J. Z. Gu, J. M. Dong, and X. Z. Zhang, Phys. Rev. C 83, 054305 (2011).
 Grasso and Anguiano (2013) M. Grasso and M. Anguiano, Phys. Rev. C 88, 054328 (2013).
 Reinhard and Nazarewicz (2010) P.G. Reinhard and W. Nazarewicz, Phys. Rev. C 81, 051303 (2010).
 Fattoyev and Piekarewicz (2011) F. J. Fattoyev and J. Piekarewicz, Phys. Rev. C 84, 064302 (2011).
 Erler et al. (2012a) J. Erler, N. Birge, M. Kortelainen, W. Nazarewicz, E. Olsen, A. Perhac, and M. Stoitsov, Nature 486, 509 (2012a).
 Gao et al. (2013) Y. Gao, J. Dobaczewski, M. Kortelainen, J. Toivanen, and D. Tarpanov, Phys. Rev. C 87, 034324 (2013).
 Kortelainen et al. (2013) M. Kortelainen, J. Erler, W. Nazarewicz, N. Birge, Y. Gao, and E. Olsen, Phys. Rev. C 88, 031305 (2013).
 Reinhard and Nazarewicz (2013) P.G. Reinhard and W. Nazarewicz, Phys. Rev. C 87, 014324 (2013).
 Reinhard et al. (2013) P.G. Reinhard, J. Piekarewicz, W. Nazarewicz, B. K. Agrawal, N. Paar, and X. RocaMaza, Phys. Rev. C 88, 034325 (2013).
 Nazarewicz et al. (2014) W. Nazarewicz, P.G. Reinhard, W. Satuła, and D. Vretenar, Eur. Phys. J. A 50, 20 (2014).
 Engel et al. (1975) Y. M. Engel, D. M. Brink, K. Goeke, S. J. Krieger, and D. Vautherin, Nucl. Phys. A 249, 215â238 (1975).
 Dobaczewski and Dudek (1996) J. Dobaczewski and J. Dudek, Acta Phys. Pol. B 27, 45 (1996).
 Perlińska et al. (2004) E. Perlińska, S. G. Rohoziński, J. Dobaczewski, and W. Nazarewicz, Phys. Rev. C 69, 014316 (2004).
 Dobaczewski et al. (2002) J. Dobaczewski, W. Nazarewicz, and M. V. Stoitsov, Eur. Phys. J. A 15, 21 (2002).
 Bertsch et al. (2009) G. F. Bertsch, C. A. Bertulani, W. Nazarewicz, N. Schunck, and M. V. Stoitsov, Phys. Rev. C 79, 034306 (2009).
 Stoitsov et al. (2003) M. V. Stoitsov, J. Dobaczewski, W. Nazarewicz, S. Pittel, and D. J. Dean, Phys. Rev. C 68, 054312 (2003).
 Schwierz et al. (2007) N. Schwierz, I. Wiedenhover, and A. Volya, arXiv:0709.3525 (2007).
 Oros (1996) A. Oros, Study of the coupling between singleparticle and collective degrees of freedom in mediummass spherical nuclei, Ph.D. thesis, University of Köln (1996).
 Alex Brown (1998) B. Alex Brown, Phys. Rev. C 58, 220 (1998).
 Chabanat et al. (1998) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
 Koopmans (1934) T. Koopmans, Physica 1, 104 (1934).
 Schunck et al. (2010) N. Schunck, J. Dobaczewski, J. McDonnell, J. Moré, W. Nazarewicz, J. Sarich, and M. V. Stoitsov, Phys. Rev. C 81, 024316 (2010).
 Satuła et al. (2012) W. Satuła, J. Dobaczewski, W. Nazarewicz, and T. Werner, Phys. Rev. C 86, 05416 (2012).
 Uozumi et al. (1994) Y. Uozumi, N. Kikuzawa, T. Sakae, M. Matoba, K. Kinoshita, S. Sajima, H. Ijiri, N. Koori, M. Nakano, and T. Maki, Phys. Rev. C 50, 263 (1994).
 Burrows (2008) T. Burrows, Nuclear Data Sheets 109, 1879 (2008).
 (65) see Supplemental Material.
 Rutz et al. (1998) K. Rutz, M. Bender, P.G. Reinhard, J. Maruhn, and W. Greiner, Nucl. Phys. A 634, 67 (1998).
 Stoitsov et al. (2013) M. V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, and S. Wild, Comput. Phys. Comm. 184, 1592 (2013).
 Furnstahl and Hackworth (1997) R. J. Furnstahl and J. C. Hackworth, Phys. Rev. C 56, 2875 (1997).
 Kortelainen et al. (2010b) M. Kortelainen, R. J. Furnstahl, W. Nazarewicz, and M. V. Stoitsov, Phys. Rev. C 82, 011304 (2010b).
 Pastore et al. (2012a) A. Pastore, D. Davesne, Y. Lallouet, M. Martini, K. Bennaceur, and J. Meyer, Phys. Rev. C 85, 054317 (2012a).
 Pastore et al. (2012b) A. Pastore, M. Martini, V. Buridon, D. Davesne, K. Bennaceur, and J. Meyer, Phys. Rev. C 86, 044308 (2012b).
 Pastore et al. (2012c) A. Pastore, K. Bennaceur, D. Davesne, and J. Meyer, Int. J. Mod. Phys. E 5, 1250041 (2012c).
 Lesinski et al. (2006) T. Lesinski, K. Bennaceur, T. Duguet, and J. Meyer, Phys. Rev. C 74, 044315 (2006).
 Pastore et al. (2013) A. Pastore, D. Davesne, K. Bennaceur, J. Meyer, and V. Hellemans, Phys. Scr. T 154, 014014 (2013).
 Hellemans et al. (2013) V. Hellemans, A. Pastore, T. Duguet, K. Bennaceur, D. Davesne, J. Meyer, M. Bender, and P.H. Heenen, Phys. Rev. C 88, 064323 (2013).
 Dobaczewski et al. (1984) J. Dobaczewski, H. Flocard, and J. Treiner, Nucl. Phys. A422 (1984).
 Brandt (1999) S. Brandt, Data Analysis: Statistical and Computational Methods for Scientists and Engineers (Springer, New York, 1999).
 Seber and Wild (1989) G. A. F. Seber and C. J. Wild, Nonlinear Regression (Wiley, 1989).
 Piekarewicz et al. (2012) J. Piekarewicz, B. K. Agrawal, G. Colò, W. Nazarewicz, N. Paar, P.G. Reinhard, X. RocaMaza, and D. Vretenar, Phys. Rev. C 85, 041302 (2012).
 Lattimer (2012) J. M. Lattimer, Annu. Rev. Nucl. Part. Sci. 62, 485 (2012).
 Erler et al. (2012b) J. Erler, N. Birge, M. Kortelainen, W. Nazarewicz, E. Olsen, A. Perhac, and M. Stoitsov, J. Phys. Conf. Ser. 402, 012030 (2012b).
 Bender et al. (2006) M. Bender, G. Bertsch, and P.H. Heenen, Phys. Rev. C 73, 034322 (2006).
 Bartel et al. (1982) J. Bartel, P. Quentin, M. Brack, C. Guet, and H.B. Håkansson, Nucl. Phys. A 386, 79 (1982).
 Tondeur (1985) F. Tondeur, Nucl. Phys. A 442, 460 (1985).
 Berger et al. (1989) J.F. Berger, M. Girod, and D. Gogny, Nucl. Phys. A 502, 85c (1989).
 Nikolov et al. (2011) N. Nikolov, N. Schunck, W. Nazarewicz, M. Bender, and J. Pei, Phys. Rev. C 83, 034305 (2011).
 Smirenkin (1993) G. Smirenkin, IAEAReport, Tech. Rep. (IAEA, 1993).
 Singh et al. (2002) B. Singh, R. Zywina, and R. Firestone, Nucl. Data Sheets 97, 241 (2002).
 Möller et al. (2009) P. Möller, A. J. Sierk, T. Ichikawa, A. Iwamoto, R. Bengtsson, H. Uhrenholt, and S. Åberg, Phys. Rev. C 79, 064304 (2009).
 Schunck et al. (2012) N. Schunck, J. Dobaczewski, J. McDonnell, W. Satuła, J. Sheikh, A. Staszczak, M. Stoitsov, and P. Toivanen, Comput. Phys. Commun. 183, 166 (2012).
 McDonnell et al. (2013) J. D. McDonnell, W. Nazarewicz, and J. A. Sheikh, Phys. Rev. C 87, 054327 (2013).
 Schunck (2013) N. Schunck, Acta Phys. Pol. B 44, 263 (2013).
 Reinhard et al. (2006) P.G. Reinhard, M. Bender, W. Nazarewicz, and T. Vertse, Phys. Rev. C 73, 014309 (2006).
 Ravenhall et al. (1983) D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys. Rev. Lett. 50, 2066 (1983).
 Pieper (2003) S. C. Pieper, Phys. Rev. Lett. 90, 252501 (2003).
 Gandolfi et al. (2011) S. Gandolfi, J. Carlson, and S. C. Pieper, Phys. Rev. Lett. 106, 012501 (2011).
 Pudliner et al. (1997) B. S. Pudliner, V. R. Pandharipande, J. Carlson, S. C. Pieper, and R. B. Wiringa, Phys. Rev. C 56, 1720 (1997).
 Pudliner et al. (1996) B. S. Pudliner, A. Smerzi, J. Carlson, V. R. Pandharipande, S. C. Pieper, and D. G. Ravenhall, Phys. Rev. Lett. 76, 2416 (1996).
 Goriely et al. (2013) S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. C 88, 024308 (2013).
 Carlsson et al. (2008) B. G. Carlsson, J. Dobaczewski, and M. Kortelainen, Phys. Rev. C 78, 044326 (2008).
 Raimondi et al. (2011) F. Raimondi, B. G. Carlsson, and J. Dobaczewski, Phys. Rev. C 83, 054311 (2011).
 Dobaczewski et al. (2012) J. Dobaczewski, K. Bennaceur, and F. Raimondi, J. Phys. G 39, 125103 (2012).