Nuclear energy density optimization: Large deformations M. Kortelainen,1,2 J. McDonnell,1,2 W. Nazarewicz,1,2,3 P.-G. Reinhard,4 J. Sarich,5 N. Schunck,6,1,2 M. V. Stoitsov,1,2 and S. M. Wild5 1Department of Physics and Astronomy, University of Tennessee, Knoxville, TN 37996, USA 2Physics Division, Oak Ridge National Laboratory, P.O. Box 2008, Oak Ridge, TN 37831, USA 3Institute of Theoretical Physics, Warsaw University, ul. Hoz˙a 69, PL-00681, Warsaw, Poland 4Institut fu¨r Theoretische Physik, Universita¨t Erlangen, D-91054 Erlangen, Germany 5Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439, USA 6Physics Division, Lawrence Livermore National Laboratory, Livermore, CA 94551, USA A new Skyrme-like energy density suitable for studies of strongly elongated nuclei has been de- 2 termined in the framework of the Hartree-Fock-Bogoliubov theory using the recently developed 1 model-based, derivative-free optimization algorithm pounders. A sensitivity analysis at the opti- 0 mal solution has revealed the importance of states at large deformations in driving the parameter- 2 ization of the functional. The good agreement with experimental data on masses and separation energies,achievedwiththepreviousparameterizationunedf0,islargelypreserved. Inaddition,the n new energy density unedf1 gives a much improved description of the fission barriers in 240Pu and a neighboring nuclei. J 0 PACSnumbers: 21.60.Jz, 21.10.-k,21.30.Fe,21.65.Mn,24.75.+i 2 ] h I. INTRODUCTION actinide regionby consideringthe experimentalinforma- t tion on the fission barrier of 240Pu. However, the op- - l timization was not performed directly at the deformed c One of the focus areas of the UNEDF SciDAC collab- HFB level; instead, a semi-classical approach was used u oration [1, 2] has been the description of the fission pro- basedontheThomas-Fermiapproximationtogetherwith n cess within a self-consistent framework based on nuclear [ shell-correction techniques. The D1S parameterization density functional theory (DFT). Until now, attempts ofthefinite-rangeGognyforcewasalsofine-tunedtothe 2 at going beyond the macroscopic-microscopic methods firstbarrierheightof240Pu[30],consideringa rotational v [3] have been carried in the context of the original self- correction to the energy of the deformed state. How- 4 consistentnuclearmean-fieldtheory[4]withSkyrme(see, ever, this fine-tuning again was not done directly at the 4 e.g., Refs. [5–10]), Gogny (see, e.g., Refs. [11–14]), and 3 HFB level but by a manual readjustment of the surface relativistic (see, e.g., Refs. [5, 15, 16]) energy density 4 coefficient of the EDF using a phenomenological model. . functionals (EDFs). The fundamental assumption of Also, in the Bsk14 EDF of the HFB-14 mass model [31] 1 the nuclear DFT is that one can describe accurately 1 by the Bruxelles-Montr´eal collaboration, data on fission a broad range of phenomena in nuclei, including ex- 1 barrierswereutilizedtooptimizetheEDFparametersby 1 cited states and large-amplitude collective motion, by adding phenomenological collective corrections, includ- : enriching the density dependence of the functional while ing a rotational one. One may, therefore, conclude that v staying at the single-reference Hartree-Fock-Bogoliubov i no EDF has ever been systematically optimized at the X (HFB) level. In this picture, beyond-mean-field correc- deformedHFB level(andwithoutphenomenologicalcor- r tions are implicitly built-in. Preliminary studies aimed rections added) by explicitly considering constraints on a atre-examiningtheoldproblemofrestoringbrokensym- states at large deformations. metries in this context are promising [17, 18]. In a previous study [25], we applied modern opti- A common challenge to both the self-consistent mean mization and statistical methods, together with high- field and the DFT approach is the need to carefully op- performance computing, to carry out EDF optimization timize the EDF parameters to the preselected pool of at the deformed HFB level, namely, the approximation observables [4, 19–25]. In particular, special attention level where the functional is later applied. The result- must be paid to optimize the parameters in the same ing EDF parameterization unedf0 yields good agree- regime where the theory will later be applied and to mentwithexperimentalmasses,radii,anddeformations. choose the fit observables accordingly. In a recent work The present work represents an extension of [25] to the [26], we showed that existing Skyrme EDFs exhibit a problem of fission. In particular, it builds on the results significant spread in bulk deformation properties, and reported in Ref. [26], which concluded that the data on re-emphasized [27, 28] that the resulting theoretical un- stronglydeformednuclearstates shouldbe consideredin certainties could be greatly reduced by considering data the optimization protocol to constrain the surface prop- corresponding to large deformation in the optimization erties of the functional. process. Let us recall that the early Skyrme-type EDF Here we propose the new EDF Skyrme parameteriza- SkM∗ [29] was in fact optimized for fission studies in the tion, unedf1, which is obtained by adding to the list 2 of fit observables the experimental excitation energies of details. (Forbrevity,wehaveomittedtheexplicitdepen- fission isomers in the actinides. To ensure that the func- dence of the densities on the coordinate r.) The J den- t tional can be used in fission and fusion studies, we have sity is the vector part of the spin-current density tensor removed the center-of-mass (c.o.m.) correction in the J . (As in our previous work [25], non-vector compo- µν spirit of the DFT. As in the case of unedf0, a sensitiv- nents of J were disregarded.) The coupling constants µν ityanalysishasbeenperformedatthesolutioninorderto are real numbers, except for Cρρ, which is taken to be t identify possible correlations between model parameters density-dependent: and assess the impact of the new class of fit observables on the resulting parameterization. Cρρ =Cρρ+Cρρ ργ. (3) t t0 tD 0 Thispaperisorganizedasfollows. InSec. IIwebriefly reviewthe theoreticalframework,establishthe notation, All volume coupling constants (Ctρρ and Ctρτ) can be and justify the removal of the c.o.m. correction. Section related to the constants characterizing the infinite nu- III defines the set of fit observables, discusses numerical clear matter [25], and this relation was used during the precisionandimplementation, andpresents the new un- optimization in order to define the range of parameter edf1 parameter set together with the results of the sen- changes. sitivity analysis. To assess the impact of fission-isomer The Coulomb contribution is treated as usual by as- data, we compare unedf1 with unedf0 in Sec. IV. In suming a point proton charge. The exchange term was Sec. V we study the performance of unedf1 with re- computed at the Slater approximation: spect to global nuclear observables, spectroscopic prop- erties, fission, and neutron droplets. Section VI contains ECoul(r)=−3e2 3 1/3ρ4/3. (4) the main conclusions and lays out future work. Exc 4 (cid:18)π(cid:19) p For the pairing energy density χ˜(r), we use the mixed II. THEORETICAL FRAMEWORK pairing description of [34] with Vq 1ρ (r) A. Time-Even Skyrme Energy Density Functional χ˜(r)= 0 1− 0 ρ˜2(r), (5) 2 (cid:20) 2 ρ (cid:21) qX=n,p c In the nuclear DFT, the total binding energy E of the where ρ˜ is the local pairing density. The value ρ =0.16 nucleus is a functional of the one-body density ρ and c fm−3 is used throughout this paper. We allow for differ- pairing ρ˜ matrices. In its quasi-local approximation, it entpairingstrengthsforprotons(Vp)andneutrons(Vn) can be written as a 3D spatial integral: 0 0 [35]. A cut-off of E = 60MeV was used to truncate cut the quasi-particle space [36]. To prevent the collapse of E[ρ,ρ˜] = d3r H(r) pairing,weusedthe Lipkin-Nogamiprocedureaccording Z to [37]. = d3r EKin(r)+χ (r)+χ (r) 0 1 Z (cid:2) +χ˜(r)+EDCiorul(r)+EECxocul(r) , (1) B. Treatment of the Center of Mass (cid:3) where H(r) is the energy density that is quasi-local (it The success of the self-consistent mean-field theory is, usually depends on derivatives with respect to the local toagreatextent,due totheconceptofsymmetrybreak- densities),time-even,scalar,isoscalar,andreal. Itisusu- ing. Aclassicexampleisthebreakingofthetranslational ally broken down into the kinetic energy (EKin(r)) and invariance by the mean field that is localized in space. nuclearpotential(forboththeparticle-holeandparticle- The associated c.o.m. correction to the binding energy particle channels, χ0,1(r) and χ˜(r), respectively) and [4,38,39],−hPˆ2 i/(2mA),isusuallyaddedtotheDFT Coulomb terms (ECoul(r) and ECoul(r)). For Skyrme c.m. Dir Exc binding energy in (1). This correction contributes typ- functionals,theparticle-holeenergydensityχ (r)+χ (r) 0 1 ically a few MeV to the total energy. Moreover, it was splits into χ (r), depending only on isoscalar densities, 0 shown that adopting approximations to this correction and χ (r), depending on isovector densities (and the 1 during the optimization of the functional could lead to isoscalarparticledensitythroughthedensitydependence significantly different surface properties [4, 17, 40]. ofthecouplingconstantCρρ;seebelow)[4,32,33]. Each 1 Since the c.o.m. correction is not additive in particle term takes the generic form number,itcausesseriousconceptualproblemswhendeal- χ (r) = Cρρρ2+Cρτρ τ +CJ2J2 ing with fission or heavy-ion fusion, that is, when one t t t t t t t t considers the split of the nucleus into several fragments, +Cρ∆ρρ ∆ρ +Cρ∇Jρ ∇·J , (2) or formationof the compound nucleus through a merger t t t t t t of two lighter ions. In fission studies, it was shown that where ρ , τ , and J (t = 0,1) can all be expressed in thecontributionofthec.o.m.correctionbetweenthetwo t t t terms of full density matrix ρ (rσ,r′σ′); see Ref. [4] for prefragmentscould amount to severalMeV near scission t 3 [41–43]. Moreover,properlycomputing this relativecon- TABLEI:Experimentalexcitation energiesoffission isomers tribution is difficult, as it reflects the degree of entan- [50] (in MeV) considered in theunedf1dataset. glement between prefragments [43, 44]. Time-dependent Hartree-Fock calculations of low-energy heavy-ion reac- Z N E tions are evenmore problematic,as there is currently no 92 144 2.750 solution to the discontinuity of the c.o.m. correction be- 92 146 2.557 tween the target+projectile system and the compound 94 146 2.800 96 146 1.900 nucleus [45, 46]: such calculations usually neglect the c.o.m. term altogether, even though EDFs employed in suchcalculationshavebeenusuallyfittedwiththec.o.m. correction included. There are, however, some excep- tions,see,e.g.,Ref.[47]. Notethatthesameproblemoc- and four excitation energies of FIs. For the FIs, we used curs with the so-calledrotationalcorrectionarising from the weight wi =0.5MeV in the χ2 objective function the breaking of rotational invariance by deformed mean fields [41, 48]. s (x)−d 2 χ2(x)= i i . (6) Another undesired property of the c.o.m. correction is (cid:18) w (cid:19) Xi i thatitslightlybreaksthevariationalnatureofHFBwhen adding or subtracting a particle [49], i.e., it violates the The χ2 weights for binding energies, proton rms radii, Koopmans theorem. The resulting s.p. energy shifts are andodd-evenmass(OEM)staggeringarethe same asin quite significant and they are of the order of the mass Ref. [25]. polarization effect related to the fact that when adding Twoassumptionsmadein[25]werealsoadoptedhere: or subtracting a particle to a closed spherical core, the (i) since the isovector effective mass cannot be reliably resulting nucleus becomes deformed. If a spherical sym- constrained by the current data, it was set to 1/M∗ = metry is imposed, the mass polarization effect is a self- 1.249 as in unedf0 and the SLy4 parameterization [v51], consistent rearrangement of all nucleons, when an odd which was the initial starting point in our optimization; particle is introduced to the system. This corresponding and (ii) since tensor terms are mostly sensitive to the energyshiftisE /A2,i.e.,about0.4MeVin40Ca[49]. Kin single-particle(s.p.) shellstructure,whichisnotdirectly As discussed in Sec. I, the EDF is supposed to cap- constrained by the unedf1 dataset, the tensor coupling ture all the physics of interest at the HFB level. In constants CJ2 and CJ2 were set to zero. In summary, other words, the functional is to be built from the full comparedwi0thunedf10[25],theoptimizationof unedf1 single-reference density matrix. While the HFB vacuum is characterized by the following: breaksthetranslationalsymmetry,theassociatedcorrec- tiontermshouldbe absorbedinthe density dependence. • The same 12 EDF parameters to be optimized, In particular, it should be possible to express the c.o.m. namely, ρ , ENM/A, KNM, aNM, LNM, M∗, Cρ∆ρ, c sym sym s 0 term not as an explicit function of A, as currently being Cρ∆ρ, Vn, Vp, Cρ∇J, and Cρ∇J; done, but through a density functional. Until a simple 1 0 0 0 1 prescriptionisproposed,however,itismoreconsistentto • 7 additional data points: 3 new masses and 4 FI simply drop all corrections that originate from a Hamil- energies with the weights w =0.5MeV; tonian view of the problem, including the c.o.m. term. And, since our focus in on fission, this is precisely what • Neglect of the c.o.m. correction term. we have done in this work. B. Numerical Precision and Implementation III. OPTIMIZATION AND SENSITIVITY ANALYSIS All HFB calculations were run with the code hfbtho [52]. The code expands the HFB solutionsonthe axially A. Experimental Dataset symmetric, deformed harmonic oscillator (HO) basis. In the optimization of unedf0, we used a spherical basis Since the focus of this study is the construction of with 20 HO shells, which was found to give a good com- an EDF optimized for fission, our experimental dataset promise between the numerical precision and computa- has been expanded by including the excitation energies tional performance. The current optimization includes (bandheads) of four fission isomers (FIs) listed in Ta- states with much larger deformation than in the ground ble I. The ground-state (g.s.) binding energies of 238U, state,andthedependence ofthe energieswithrespectto 240Pu, and 242Cm, not included in the previous unedf0 the set of basis states is more significant. dataset [25], were added for consistency (the g.s. energy Intheunedf1optimization,allquantitiesbutthefour of236Uwasalreadythere). Consequently,comparedwith fissionisomerswerecomputedwiththesphericalHOba- unedf0, the unedf1 dataset contains seven new data sis of N = 20 shells, which includes N = 1771 basis sh points: three additional g.s. masses of deformed nuclei, states. For the fission isomers, we adopted a stretched 4 HO basis with deformation β = 0.4. The basis con- hfbtho runs over a more traditional, derivative-free tains up to N = 50 oscillator shells with an upper Nelder-Mead optimization method [53]. Hence, with- sh limit of N = 1771 basis states with lowest HO s.p. en- outthealgorithmicandcomputationaladvancementsde- ergies. The oscillator frequency ω3 = ω2ω was set at tailedabove,asimilaroptimizationcouldhavepreviously 0 ⊥ k ~ω = 1.2×41/A1/3MeV. As seen in Fig. 1, at this se- consumedamonthofcomputationsusing80coresofthe 0 lection of the HO basis, the dependence of FI energies Fusion cluster. on the basis deformation remains fairly constant around We emphasize that, strictly speaking, both the un- β =0.4. Moreover,therangeofvariationsissignificantly edf0andtheunedf1parameterizationsobtainedinthis less than the corresponding χ2 weight, w =0.5MeV. workshouldalwaysbeusedintheiroriginalenvironment. i In particular, the pairing EDF should be that of Eq. (5) usedwiththeoriginalpairingspacecutoff;pairingcalcu- $,-./0 !"#$ !&’() lationsmustbecomplementedbytheLipkin-Nogamipre- scription; and the proton and neutron pairing strengths C B !"%$ !&!*+ must not vary from the values determined by our opti- 6 A mization. In short, contrary to usual practice, there is 5@ no flexibility in the treatment of the pairing channel. ? > 9 6 ; - C. Result of the Optimization: UNEDF1 5 ; Parameter Set 8 :4 2 =4: The startingpoint forour poundersoptimizationwas < thepreviouslyobtainedunedf0parameterization. After - 177 simulations, the algorithm reached the new optimal result. The resulting parameter set is listed in Table II. Thefirstsixparameterswererestrictedtoliewithinfinite 123435.6789+2:48;5 !5 bounds, also listed in Table II, that were not allowed to beviolatedduringtheoptimizationprocedure. Ascanbe FIG. 1: (Color online) Excitation energies of fission isomers seen,parametersENM/AandKNM areonthe boundary consideredintheunedf1optimizationasfunctionsoftheHO value. In the case of unedf0, we recall that KNM and basis deformation. 1/M∗ also ended up at their respective boundaries. The s saturation density ρ is given with more digits than the c Optimization calculations were performed on Ar- other parameters. Such extra precision is needed when gonne National Laboratory’s Fusion cluster, managed computing volume coupling constants [25]. by Argonne’s Laboratory Computing Resource Center (LCRC). Fusion consists of 320 computing nodes, each TABLE II: Optimized parameter set unedf1. Listed are with dual quad-core Pentium Xeon processors. By us- bounds used in the optimization, final optimized parameter ing Intel’s Math Kernel Library and the Intel Fortran values, standard deviations, and 95% confidenceintervals. compiler (ifort), we were able to run hfbtho in almost half the time when comparedwith the prebuilt reference x Bounds xˆ(fin.) σ 95% CI BLAS library implementation and GNU’s gfortran com- ρ [0.15,0.17] 0.15871 0.00042 [ 0.158, 0.159] c piler. We werealsoable to dramaticallyreduce the wall- ENM/A [-16.2,-15.8] -15.800 – – clocktimeofanhfbthocomputationbyusingOpenMP KNM [220, 260] 220.000 – – atthe nodeleveltoparallelizekeycomputationalbottle- aNsyMm [28, 36] 28.987 0.604 [ 28.152, 29.822] necks. These bottlenecks involved iteratively computing LNsyMm [40, 100] 40.005 13.136 [ 21.841, 58.168] the eigenvalues and eigenvectors of the (Ω, π) blocks of 1/Ms∗ [0.9, 1.5] 0.992 0.123 [ 0.823, 1.162] the HFB matrix, as well as the density calculations re- C0ρ∆ρ [−∞,+∞] -45.135 5.361 [ -52.548, -37.722] flecting the same block pattern. OpenMP allowed us to C1ρ∆ρ [−∞,+∞] -145.382 52.169 [-217.515, -73.250] Vn [−∞,+∞] -186.065 18.516 [-211.666,-160.464] dynamicallyassignprocessorstoblocksofdataforparal- 0 Vp [−∞,+∞] -206.580 13.049 [-224.622,-188.538] lel processing, which further reduced the wall-clock time 0 Cρ∇J [−∞,+∞] -74.026 5.048 [ -81.006, -67.046] by a factor of 6 when running on an eight-core node. 0 Cρ∇J [−∞,+∞] -35.658 23.147 [ -67.663, -3.654] The parameter estimation computations presented in 1 this paper ran 218 total simulations of hfbtho for each nucleus in the dataset, using 80 compute nodes (640 cores) for 5.67 hours. As highlighted in [25], using the We first note that the same minimum was obtained pounders algorithm (Practical Optimization Using No by starting either from the unedf0 solution or from the Derivatives (for Squares)) on the type of fitting prob- unedf1ex parameterization discussed below: this gives lem considered here requires more than 10 times fewer us confidence that the parameter set listed in Table II 5 is sufficiently robust. We can then observe that most of ∗ M s UNEDF1 the parameter values of unedf1 are fairly close to those 1/ of unedf0 [25]. There are, nevertheless, a couple of no- ρ∇J 1 C1 table exceptions. First, the magnitude of Cρ∆ρ is now ρ∇J 1 C0 much larger. This is potentially dangerous, as it might p 0.75 V0 triggerscalar-isovectorinstabilitiesinthefunctionalthat n could appear in neutron-rich nuclei [54, 55]. (Our mass- V0 table calculations with unedf1 do not show indications ρ∆ρ 0.5 of instability in even-even nuclei.) Second, C1ρ∇J has Cρ∆1ρ drifted considerably from its initial value, even changing C0 NM m 0.25 sign. These two coupling constants control the isovector Lsy M surfacepropertiesofthenucleus;hence,onlypropercon- N m asy straintsontheshellstructurelike,forexample,spin-orbit 0 ρc splitting in neutron-rich nuclei will allow these terms to bsteanptisnnaeredrdeolawtinv.elyFourncthonesmtroamineendt,,absoetvhidceonucpedlinaglsocobny- ρc aNsMym LNsMym Cρ∆0ρ Cρ∆1ρ Vn0 Vp0 Cρ∇0J Cρ∇1J1/M∗s their relatively large σ-value shown in Table II. FIG.2: (Coloronline) Absolutevaluesofthecorrelation ma- D. Sensitivity Analysis trix of Table III presented in a color-coded format. 1. Correlation Matrix of unedf1 each data type. Let us recall that the sensitivity matrix S is defined as We have performed a sensitivity analysis at the solu- S(x)= J(x)JT(x) −1J(x), (7) tion of the optimization. All residual derivatives were estimated by using the optimal finite-difference proce- where J(x) is the Jaco(cid:2)bianmatrix(cid:3). The results areillus- dure detailed in Ref. [56]. Since some of the parameters trated in Fig. 3, where we have summed absolute values ranattheirboundsduringtheoptimization,thesensitiv- ofeachdatatypeforeachparameter. Thetotalstrengths ity analysis was carried out in a subspace that does not for each parameter were then normalized to 100%. Note contain these parameters. The same strategy was also that the fission isomer excitation energies represent less used in the previous sensitivity analysis of the unedf0 than 4% of the total number of data points but account parameterization; we refer to [25] for a detailed discus- fortypically30%ofthevariationoftheparameterset. In sion of the available options in the case of constrained thecaseofthesymmetryenergycoefficient,thispercent- optimization. In Table II we list the standard deviation, age is even 75% (see Sec. IV for more discussion). Com- σ, and 95% confidence interval (CI) for each parameter paredwith unedf0, we find that the overalldependence at the solution. As discussed in Sec. IIIC, the standard ontheprotonradiihassignificantlydecreased,exceptfor deviations of most of the parameters is relatively small. ρ and Cρ∆ρ, and that the dependence on the OES has The correlation matrix for the unedf1 parameter set actually1increased. This kind of analysis, however, does is presented in Table III. It was calculated as in [25] not address the importance of a particular data point to and, similarly, corresponds to the 10-dimensional sub- the obtained optimal solution. space of the parameters that are not at their boundary A complementary way to study the impact of an indi- value. Generally,mostoftheparametersareonlyslightly vidual datum on the obtained parameter set is therefore correlated to each other, with a few notable exceptions presented in Fig. 4. Here, we have plotted the amount (correlationsbelow0.8arenotverysignificantfromasta- of variation tisticalviewpoint). Thestrongcorrelationbetween1/M∗ s andbothVn andVp hadalreadybeennoticedintheun- δx 2 edfnb para0meter s0et of Ref. [25] and reflects the inter- ||δx/σ||=vu (cid:18) σk(cid:19) (8) play between the level density at the Fermi-surface and utXk the size of pairing correlations. Similarly, both pairing for the optimal parameter set when data points d are i strengths are strongly correlated with Cρ∆ρ, which can changed by an amount of 0.1w one by one. As can be 0 i alsoberelatedtosurfacepropertiesofthefunctional. In- seen,thevariationsaresmalloverall,assuringusthatthe terestingly, both pairing strengths are now strongly cor- dataset was chosen correctly. The masses of the double related with one another, which was not the case with magicnuclei208Pband58Niseemtohavethebiggestrel- unedf0. The same correlation matrix of Table III is ativeimpactonthe optimalparameterset. Onecanalso shown graphically in Fig. 2. see that the sensitivity of the parameters on the new FI Next, westudy the overallimpactofeachdatatype in dataislargerthantheaveragedatumpoint. Bycontrast, ourχ2 functionontheobtainedparameterset. Asin[25] the dependence ofthe parameterizationonthe massesof wecalculatethepartialsumsofthesensitivitymatrixfor deformed actinide and rare earth nuclei is weaker. 6 TABLE III: Correlation matrix for unedf1parameter set. Correlations greater than 0.8 (in absolute value) are in boldface. ρ 1.00 c aNM -0.35 1.00 sym LNM -0.14 0.71 1.00 sym 1/M∗ 0.32 0.23 0.36 1.00 s Cρ∆ρ -0.25 -0.25 -0.35 -0.99 1.00 0 Cρ∆ρ -0.06 -0.15 -0.77 -0.22 0.19 1.00 1 Vn -0.32 -0.22 -0.36 -0.99 0.98 0.22 1.00 0 Vp -0.33 -0.18 -0.29 -0.97 0.97 0.15 0.96 1.00 0 Cρ∇J -0.14 -0.20 -0.32 -0.86 0.91 0.22 0.85 0.84 1.00 0 Cρ∇J 0.05 -0.17 -0.13 -0.10 0.07 0.21 0.10 0.07 -0.03 1.00 1 ρ aNM LNM 1/M∗ Cρ∆ρ Cρ∆ρ Vn Vp Cρ∇J Cρ∇J c sym sym s 0 1 0 0 0 1 Masses (sph) Charge radii OESFI Masses (def) )*+’ !"#$%$&’( ,*+#-( ./ 0% 2$1* FIG.3: (Coloronline)Sensitivityofunedf1todifferenttypes FIG. 4: (Color online) Overall change in x for the unedf1 of data entering theχ2 function. when thedatumdi is changed byan amount of 0.1wi oneby one. ThefourrightmostdatapointsmarkedFIcorrespondto excitation energies offission isomers. Theresultsfor unedf0 and unedf1ex of Sec. IIID2 are also shown. 2. Discussion of the Coulomb Exchange Term It has been argued in [57] that removing the Coulomb this additional parameter α is given in Table IV. Our ex exchange term from the functional could improve the objective function is slightly decreased from 51.058 to overall fit on nuclear binding energies. This procedure 49.341 when this term is present. Overall, both param- hadbeenmotivatedbytheearlierworksRefs.[58,59],in eterizations, with and without the Coulomb exchange which similar ameliorations, albeit on a smaller dataset, screening term, are very similar. However, one can see wereobserved. The originofthese adhoc manipulations that the 95% CI is relatively large for the screening was the observation that many-body effects induced by parameter, the value of which is also close to 1 (full the long-range Coulomb force among protons manifest Coulomb exchange). We recall that this confidence in- themselves in the form ofa (positive) correlationenergy, tervalisextractedfromthe correlationmatrixcomputed which, to some extent, can cancel out the (negative) ex- in the 10-D space of “inactive” parameters, namely, the changeterm[60,61]. Since suchanexchange-correlation space of the 10 parameters that are not at their bound effectisabsentfromthestandardSkyrmefunctional,one andthus activelyconstrained. Ifonecomputes the Jaco- could feel justified to simulate it by effectively screen- bian matrix in the original 13-D space of all parameters ing the Coulombexchangetermwith anempiricalfactor with a tangent plane approximation to account for the 0≤αex ≤1. Thespecialcasesαex =0andαex =1give, three active parameters ENM, KNM, and LNM, we find sym respectively, the case without and with full Coulombex- that the 95% CI for the screening parameter becomes change. [-1.663,3.290]. This implies that α is basically notcon- ex The result of the optimization of the functional with strained with the current dataset. 7 especially if the coupling constants driving shell effects TABLE IV: Optimized parameter set unedf1ex. Listed are are somewhat constrained by the dataset. In the lan- boundsusedintheoptimization, finaloptimized parameters, guage of the leptodermous expansion of Sec. IVC, this standard deviations, and 95% confidenceintervals. implies that both the surface and surface-symmetry en- x Bounds xˆ(fin.) σ 95% CI ergycoefficients(whichdependinanontrivialwayonthe ρc [0.15, 0.17] 0.15837 0.00049 [ 0.158, 0.159] couplingconstantsofthefunctional)shouldbeimpacted. ENM/A [-16.2, -15.8] -15.800 – – On the other hand, we may assume that the isospin de- KNM [220, 260] 220.000 – – pendence of the binding energy (i.e., the total symmetry aNsyMm [28, 36] 28.384 0.711 [ 27.417, 29.351] energy) is relatively well constrained by the several long LNsyMm [40, 100] 40.000 – – isotopic sequences present in our dataset. We therefore 1/M∗ [0.9, 1.5] 1.002 0.123 [ 0.835, 1.169] s seethatthe requirementofhavingthe fullsymmetryen- Cρ∆ρ [−∞,+∞] -44.602 5.349 [ -51.872, -37.331] 0 ergy constrained together with a relatively large varia- Cρ∆ρ [−∞,+∞] -180.956 47.890 [-246.050,-115.863] 1 tionofthe surfaceterms shouldleadto a relativelylarge Vn [−∞,+∞] -187.469 18.525 [-212.649,-162.288] 0 variation of the volume symmetry aNM, which is indeed Vp [−∞,+∞] -207.209 13.106 [-225.024,-189.395] sym C0ρ∇J [−∞,+∞] -74.339 5.187 [ -81.389, -67.289] observed in Fig. 4. C0ρ∇J [−∞,+∞] -38.837 23.435 [ -70.690, -6.984] One can now understand the difference of behavior of 1 aNM under a change of data when the Coulomb screen- α [ 0, 1] 0.813 0.154 [ 0.604, 1.023] sym ex ingtermispresent: accordingto[60,61],themany-body Coulombcorrelationenergythatissimulatedbyα <1 ex essentially represents a proton surface effect. Changes in bulk surface properties triggered by variations in the The dependence of every parameter on the four types excitation energy of fission isomers can be entirely ab- of data included in the dataset (masses, charge radii, sorbedby areadjustmentofα ,especially sincethe lat- ex OES, and FI data) reveal an interesting consequence of ter is poorly constrained by the other data, rather than the screening of the Coulomb exchange term. Figure 5 by aNM. Lastly, we note that the Coulomb exchange sym shows the analogue of Fig. 3 when the screening pa- term, which is approximated by the usual local Slater rameter αex is included. Note the striking difference in expression, may get worse at large deformations [62]. the bar plot for the symmetry energy parameter aNsyMm: In summary, considering that (i) αex is poorly con- fluctuations in this parameter under a variation of the strained by the data, yet may affect significantly other excitation energy of the fission isomers are reduced to parameters like aNM and (ii) α does not significantly sym ex less than 10%, compared with nearly 75% when the full improve the quality of the fit, we decided to retain the Coulomb exchange term is computed. full Coulomb exchange term in the present unedf1 pa- rameterization. IV. CHARACTERIZATION OF UNEDF1 PARAMETERIZATION In this section, we discuss general properties of the unedf1 parameterizationand compare it with unedf0. A. Energy Density in (t,x) Parameterization For practical applications, it is useful to express the coupling constants of unedf0 and unedf1 in the tra- ditional (t,x)-parameterization of the standard Skyrme force, see Appendix A of [4]. The results are given in Table V. As can be seen, in the (t,x)-parameterization the two functionals are quite different. This is to be FIG. 5: (Color online) Similar as in Fig. 3 but for the un- expectedastherelationbetweentheC and(t,x)param- edf1ex parameter set. eterizations is partially nonlinear [63]. This behavior can be qualitatively understood by re- calling a few simple facts about the bulk nuclear energy B. Energy Density Parameters in Natural Units and deformation. A variationof the excitation energy of very deformed states such as fission isomers essentially TheEDFparameterscanalsobeexpressedintermsof affects the bulk surface properties of the functional— natural units [63]. In Table VI we list the parameter set 8 dropmodel (LDM), one needs to carryout the leptoder- TABLE V: Parameters (t,x) of unedf0and unedf1. mousexpansion. Thegeneralstrategybehindtheexpan- Par. unedf0 unedf1 Units sionofnuclearEDFhasbeendiscussedinRef.[64],where t −1883.68781034 −2078.32802326 MeV·fm3 one can find the relevant literature on this topic. The 0 t 277.50021224 239.40081204 MeV·fm5 starting point is the LDM binding energy per nucleon 1 t 608.43090559 1575.11954190 MeV·fm5 expanded in the inverse radius (∝ A−1/3) and neutron 2 t3 13901.94834463 14263.64624708 MeV·fm3+3γ excess I =(N −Z)/A: x 0.00974375 0.05375692 - 0 x −1.77784395 −5.07723238 - E(A,I) = a + a A−1/3 + a A−2/3 1 vol surf curv x −1.67699035 −1.36650561 - x23 −0.38079041 −0.16249117 - + asymI2 + assymA−1/3I2 (9) b 125.16100000 38.36807206 MeV·fm5 b4′ −91.2604000 71.31652223 MeV·fm5 + a(sy2)mI4. 4 γ 0.32195599 0.27001801 - For any functional, our approach consists of combin- ing nuclear matter (NM) calculations with Hartree-Fock (HF)calculationsforalargesetofsphericalnucleitoex- tract by linear regression the various parameters of the of unedf1 in standard units and in natural units. Here expansion (9) according to the following procedure. we have used the same value for the scale Λ=687MeV, characterizing the breakdown of the chiral effective the- ory,whichhasbeenfoundinRef.[63]. Fromthenumbers TABLEVII:Liquiddropcoefficients of unedf0and unedf1 inTableVI onecanseethatmostoftheunedf1 param- (all in MeV). eters are natural, with only two minor exceptions. First, a a a(2) a a a becausetheeffectivemassMs∗ inouroptimumiscloseto unedf0 -16v.0ol56 30s.y5m43 4.s4y1m8 1s8u.r7f 7cu.1rv -s4sy4m unity, the Cρτ is abnormally small. Second, Cρ∆ρ seems unedf1 -15.800 28.987 3.637 16.7 8.8 -29 0 1 beontheborderlineofbeingunnaturallylarge. AsTable II indicates, however,the standarddeviationfor this pa- First, the bulk parameters a and a are directly rameterisratherlarge. Ithastobenoted,however,that vol sym obtained from NM calculations. Second, the smooth en- there is nothing unusual about the magnitude of C1ρ∆ρ. ergy per nucleon E¯(A,I) is extracted from the spherical Indeed, some examples of EDF parameterizations with HF calculations of (A,I) nuclei by removing the shell similarorlargervaluesofCρ∆ρ canbe found inFig.2of 1 correction [64]. The isoscalar coefficients of the expan- [63]. sion(9) canthen be deduced fromthe smoothenergyby plotting TABLE VI: Coupling constants of unedf0 and unedf1 in E¯(A,0)−a A1/3 −→a +a A−1/3 (10) vol surf curv normal units and in natural units. Value Λ = 687MeV was (cid:2) (cid:3) used. asafunctionofA−1/3. Thea coefficientisobtainedas surf unedf0 unedf1 the extrapolation of the curve to A−1/3 −→0. The cur- Coupling Normal Natural Normal Natural vature coefficient a is then estimated from the slope curv Constant Units Units Units Units of the line. C0ρ0ρ -706.38 -0.795 -779.37 -0.878 The determination of isovector coefficients starts with CC1ρρ0ρρ 284608..2867 00..297011 288981..0418 00..392347 the second-order symmetry coefficient a(sy2)m. It is easily 0D estimatedbysystematiccalculationsinasymmetricNM. Cρρ -69.77 -0.072 -201.37 -0.212 1D Defining Cρτ -12.92 -0.176 -0.99 -0.014 0 C1ρτ -45.08 -0.616 -33.52 -0.458 a(eff)(∞,I) = ENM(∞,I)−ENM(∞,0) /I2 Cρ∆ρ -55.26 -0.755 -45.14 -0.616 sym (11) C0ρ∆ρ -55.62 -0.759 -145.38 -1.985 −→ (cid:2)asym+a(sy2)mI2, (cid:3) 1 C0ρ∇J -79.53 -1.086 -74.03 -1.011 one can extract the second-order symmetry coefficient C1ρ∇J 45.63 0.623 -35.66 -0.487 from the slope of a(eff)(∞,I) vs. I2. Extracting the γ 0.3220 0.2700 sym surface-symmetry coefficient is more involved. We first introduce the effective symmetry coefficient for a finite nucleus as a(eff)(A,I) = E¯(A,I)−E¯(A,0) /I2 sym (12) C. Leptodermous Expansion −→ (cid:2)a +a A−1/3(cid:3)+a(2) I2. sym ssym sym To extract global properties of the energy functional In nuclear matter (A → +∞), the effective symmetry and relate them to the familiar constants of the liquid coefficient reduces to (11). To avoid multidimensional 9 !"#$%& !"#$%’ EDFwasoptimizedprimarilytonuclearmasses. Indeed, 0 the main difference between unedf0 and unedf1 lies / -28 . -42 in surface properties. Relatively low values of asurf and - -30 a of unedf1 reflect the new constraints on the FI -44 ssym , + data and the neglect of the c.o.m. term. Again, com- 12 -46 -32 paring the LDM values of unedf1 with those in Ta- ’ ’ & ble I of Ref. [64], we note that the LDM parameters of unedf1 are closest to those of the BSk6 EDF [65] 0 17.4 /19.4 (a = 17.3MeV and a = −33MeV) and the LSD . surf ssym - 17.2 LDM [66] (asurf = 17.0MeV and assym = −38.9MeV). ,19.0 In Sec. VC, we shall see that the reduced surface en- +()* 16.8 ergy of unedf1 with respect to unedf0 has profound ’ & consequencesfor the descriptionof fissionbarriersin the 0 0.04 0.08 0 0.04 0.08 actinides. To see this reduction more clearly, we inspect !"#$% the effective surface coefficient FIG.6: (Coloronline)Surface-symmetrycoefficient(15)(up- a(eff) =a +a I2. (16) surf surf ssym perpanels)andsurfacecoefficient(lowerpanels)versusA−1/3 for unedf0(left) and unedf1(right). For 240Pu, the value of a(eff) is 15.33MeV for unedf1, surf 16.63MeV for unedf0, 15.87MeV for SLy4, 15.75MeV for BSk6, 15.15MeV for SkM∗, and 15.17MeV for LSD. regression analysis, we introduce the reduced symmetry coefficient by subtracting the I2-dependent part of the (eff) NM limit to asym: V. PERFORMANCE OF UNEDF1 a(red)(A,I)= E¯(A,I)−E¯(A,0) /I2−a(2) I2. (13) sym sym A. Global Mass Table (cid:2) (cid:3) At the perfect LDM limit, the quantity a(red)(A,I) sym One of the key elements required from the universal should not depend on the neutron excess. At small EDF is the ability to predict global nuclear properties, isospins, however, numerical uncertainties in the shell- such as masses, radii, and deformations, across the nu- correction procedure are amplified by the I2 denomina- clear chart, from drip line to drip line. We have cal- tor. Inpractice,itismoreefficienttobuildanI-averaged culated the g.s. mass table with unedf1 for even-even reduced asymmetry coefficient, nuclei with N,Z > 8. Table VIII contains the rms de- viations from experiment for binding energies, separa- 1 b a(red) = dIa(red)(A,I), (14) tion energies,averagedthree-pointodd-evenmass differ- sym b−aZa sym ences, and proton radii. Since the set of fit observables constraining unedf1 is biased toward heavy nuclei, we where we choose a = 0.1 and b = 0.2 [64]. The surface- also show rms deviations for light (A < 80) and heavy symmetry energy is then obtained from (A≥80) subsets. Figure 7 displays the binding energy residuals (i.e., a(syremd)−asym A1/3 =assym. (15) deviations from experiment). From this figure and Ta- h i ble VIII, we can see a couple of trends. First, the en- Figure 6 illustrates the numerical accuracy of the ergyresidualswithunedf1arelargerthanthoseforun- method of evaluation for the surface and surface- edf0. Thisresultisnotsurprising,asthenewdataonFI symmetrycoefficientsoftheLDM.Thedashedbluelines andtheremovalofthecenter-of-masscorrectionstrongly indicate the fitting lines from which the final values of disfavorsthe lightestnuclei during the optimizationpro- a , a , and a are deduced. The case of unedf1 cess. Second,thecharacteristicarc-likebehaviorbetween surf curv ssym seems to be clear. The trend of the surface energy for the magic numbers is pronounced, althoughthis trendis unedf0 is less clean. The two groups of nuclei, huge much weaker than, for example, for the SLy4 functional and large, seem to follow slightly different slopes, and (see Fig. 7 of Ref. [25]). the fit represents a compromise. The resulting surface In Fig. 8 we display the residuals of two-neutron and and curvature energy have to be taken with care. two-proton separation energies. Again, the emphasis of The LDM parameters of unedf0 and unedf1 are unedf1 on heavy nuclei is clearly seen, and the corre- given in Table VII. As seen in Table I of Ref. [64] and sponding rms deviations in Table VIII quantify this fea- Fig. 1 of Ref. [26], symmetry coefficients of phenomeno- ture. Noticethattwo-protonseparationenergiesaresys- logicalLDMmass models cluster arounda =30MeV tematically overestimated. The same trend is observed sym and a = −45MeV, and the unedf0 values are right for the unedf0 functional. We can speculate about ssym in the middle. This result is not surprising, as this sources for this effect: (i) Following the arguments of 10 TABLE VIII: Root-mean-square deviations from the exper- imental values for unedf0 and unedf1 for different observ- 4 UNEDF1 ablescalculatedineven-evensystems: bindingenergyE,two- ) neutronseparation energyS2n,two-protonseparationenergy V S2p, three-point odd-even mass difference ∆˜(n3) (all in MeV), Me 2 RanMdSrmdesvpiartoitoonnforraduiineRdpf0(inanfdmu).neCdofl1u,mannsdatrheeonbusmerbvaebrloef, (p 0 x e data points. E Observable unedf0 unedf1 # h-2 E 1.428 1.912 555 Et E (A<80) 2.092 2.566 113 -4 E (A≥80) 1.200 1.705 442 S 0.758 0.752 500 2n 0 20 40 60 80 100 120 140 S (A<80) 1.447 1.161 99 2n Neutron Number S (A≥80) 0.446 0.609 401 2n S 0.862 0.791 477 2p FIG. 7: (Color online) Binding energy residuals between un- S (A<80) 1.496 1.264 96 S2p (A≥80) 0.605 0.618 381 edf1 results and experiment for 555 even-even nuclei. Iso- 2p topic chains of a given element are connected by lines. ∆˜(3) 0.355 0.358 442 n ∆˜(3) (A<80) 0.401 0.388 89 n ∆˜(3) (A≥80) 0.342 0.350 353 B. Spherical Shell Structure n ∆˜(3) 0.258 0.261 395 p ∆˜(3) (A<80) 0.346 0.304 83 Thenuclearshellstructurehasasubstantialimpacton p ∆˜(3) (A≥80) 0.229 0.248 312 manynuclearproperties. Notably,thesingle-particlelev- p els close to the Fermi surface affect many nuclear prop- Rp 0.017 0.017 49 erties such as the strength of pairing correlations and R (A<80) 0.022 0.019 16 p deformability. Compared with our previous work [25], R (A≥80) 0.013 0.015 33 p the s.p. energies that we report here have been obtained from proper blocking calculations at the HFB+LN level [69], instead of being the eigenvalues of the HF Hamil- tonian. This choice is motivated by the need to stay within a logically consistent framework: both unedf0 andunedf1 havebeen optimizedatthe HFB+LN level, andhenceshouldbeemployedexclusivelyinthiscontext. [60], one may argue that the standard Skyrme function- Moreover,inthe nuclearmean-fieldtheorywitheffective als, such as unedf0 and unedf1, lack the capability to interactions,HF eigenvaluesarea poor representationof describe many-body Coulomb effects; (ii) The explicit s.p.energies;see[70,71]forarecentstudy. InaDFTap- contribution of the Coulomb field to the pairing channel proach,however,it is assumedthat the generalizedform [67, 68] is not taken into account. It is expected that of the energy density may effectively account for beyond separate pairing strengths for neutrons and protons, as mean-fieldeffectssuchasparticle-vibrationcouplings. In in unedf1, will partly account for this missing contri- addition to this theoretical argument, let us recall that bution [35]. However,the existence of nontrivialcorrela- s.p. energies are not observables but model-dependent tions between pairing strengths and other parameters of quantities extracted experimentally from binding ener- the functional (see Table III) may have consequences for gies of excited states in odd nuclei. Systematic errors observables such as two-proton separation energies. canthus be reducedbyworkingexclusivelywith binding energies. To compare unedf0 and unedf1 quantitatively, we To this end, we computed a number of one-quasi- canassesstheirperformanceonvariousobservableslisted particle(q.p.) configurationsfor the odd-massneighbors inTableVIII. Itisexpectedthatsincenewconstraintson of 16O, 40Ca, 48Ca, 56Ni, 132Sn, and 208Pb. Calcula- fissionisomershavebeenaddedwhenoptimizingunedf1 tionsweredoneattheequalfillingapproximation,which whilekeepingthesamenumberofparametersoptimized, is an excellent approximation to the full time-reversal, the rms deviations for masses and separation energies symmetry-breaking blocking scheme [69]. Blocking q.p. must increase. Indeed, the rms deviation for the masses states induces a small shape polarization [49], which in is slightly worse for unedf1, for both light and heavy turn leads to a fragmentationofsphericals.p. orbitalsof nuclei. Interestingly, the quality of S values remains angular momentum j into 2j+1 levels Ω = −j,...,+j. 2n roughly the same in both cases, as is true also for odd- In principle, the “experimental” s.p. energy should be even mass differences and proton radii. theaverageenergyoverallthe 2j+1blockingconfigura-