Nuclear energy density optimization: Large deformations

0.47 MB·English

by  M. Kortelainen

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-