A Bayesian Analysis of Nuclear Deformation Properties with Skyrme Energy Functionals
LLLNL-JRNL-810569
A Bayesian Analysis of Nuclear DeformationProperties with Skyrme Energy Functionals
N. Schunck , K. R. Quinlan , J. Bernstein Nuclear and Chemical Sciences Division, Lawrence Livermore National Laboratory,Livermore, CA 94551, USA Applied Statistics Group, Lawrence Livermore National Laboratory, Livermore, CA94551, USAE-mail: [email protected]
May 2020
Abstract.
In spite of numerous scientific and practical applications, there is stillno comprehensive theoretical description of the nuclear fission process based solely onprotons, neutrons and their interactions. The most advanced simulations of fissionare currently carried out within nuclear density functional theory ( dft ). In spite ofbeing fully quantum-mechanical and rooted in the theory of nuclear forces, dft stilldepends on a dozen or so parameters characterizing the energy functional. Calibratingthese parameters on experimental data results in uncertainties that must be quantifiedfor applications. This task is very challenging because of the high computational costof dft calculations for fission. In this paper, we use Gaussian processes to buildemulators of dft models in order to quantify and propagate statistical uncertaintiesof theoretical predictions for a range of nuclear deformations relevant to describing thefission process.PACS numbers: 21.10.-k, 21.30.Fe, 21.60.Jz, 21.65.Mn
Submitted to:
J. Phys. G: Nucl. Phys.
1. Introduction
Nuclear fission plays a key role in a number of basic and applied science problems, fromunderstanding the origin of elements in the universe [1] to the stability of superheavyelements [2]. From a fundamental science perspective, one would like to be able todescribe nuclear fission in a fully quantum-mechanical way as emerging from nuclearforces within the nucleus. Such a “microscopic” picture poses formidable challenges totheorists as it would require solving the quantum many-body problem of fermions ininteraction even though nuclear forces remain poorly known. Currently, nuclear densityfunctional theory represents our best attempt to tackle such a problem [3].1 a r X i v : . [ nu c l - t h ] J un LNL-JRNL-810569Density functional theory ( dft ) is a general approach to the quantum many-bodyproblem that is designed to scale well with particle number [4, 5]. In the particular caseof nuclear physics, effective, in-medium nuclear forces are encoded in an energy densityfunctional ( edf ), which is a function of the intrinsic density of nucleons. Given this edf ,various theoretical techniques allow computing a number of nuclear properties rangingfrom ground-state properties, low-lying excited spectra or large-amplitude collectivemotion such as fission or nuclear reactions [6].While dft has been successfully applied in many areas of nuclear science, it shouldbe viewed as an imperfect, phenomenological model. In particular, the parameters ofthe edf are unknown and must be calibrated on a set of experimental data, and thisprocess depends on the particular level of approximation within dft [7, 8, 9]. Overthe past decade, there have been numerous attempts to quantify the uncertainties ofthis calibration on predictions, but these earlier studies have focused mostly on variousground-state properties such as masses [10, 11, 12, 13], drip-line properties [14, 15], orneutron skin [16, 17]. There are still relatively few examples where either covarianceor Bayesian techniques were applied to more complex problems such as collectiveexcitations [18] or fission barriers [19]. Yet, in light of the computational cost of fissioncalculations, such analyses are essential to better identify model weaknesses.The main goal of this paper is to assess whether standard statistical methods such asBayesian inference with Gaussian processes can be used with confidence in the emulationand uncertainty quantification of fission models. Specifically, we wish to extend thework of [19] to (i) build an emulator of fission pathways and (ii) build an emulatorof the location and characteristics of scission configurations (the point where the twofragment are formed) and (iii) determine the posterior distribution of edf parametersconditioned on fission barriers .Our paper is organized in three main sections. Section 2 recalls some basic elementsof nuclear energy density functional theory for the particular case of Skyrme functionals.In section 3, we summarize the statistical models that we used in the analysis includingGaussian processes and Bayesian inference with Markov-Chain Monte-Carlo sampling.Finally, we show in section 4 a selection of results for the benchmark case of the potentialenergy curve of the
Pu nucleus.
2. Physics Background
In this section, we give a brief summary of the (single-reference) energy densityfunctional theory with Skyrme generators that will be used throughout this paper.
The general framework for this work is the Hartree-Fock-Bogoliubov ( hfb ) theory, wherethe nuclear many-body wave function takes the form of a quasiparticle vacuum [6, 20, 21].In the hfb approximation, all degrees of freedom for the system are encoded in the one-2LNL-JRNL-810569body density matrix ρ and pairing tensor κ (and its complex conjugate κ ∗ ). In particular,the total energy of the nucleus reads E [ ρ, κ, κ ∗ ] = E nuc [ ρ ] + E Cou [ ρ ] + E pair [ ρ, κ, κ ∗ ] . (1)Here, we use a Skyrme-like energy density functional ( edf ) for the nuclear part, E nuc [ ρ ] = (cid:88) t =0 , (cid:90) d r χ t ( r ) , (2)where the Skyrme edf includes the kinetic energy term and reads χ t ( r ) = C ρρt ρ t + C ρτt ρ t τ t + C JJt J t + C ρ ∆ ρt ρ t ∆ ρ t + C ρ ∇ Jt ρ t ∇ · J t . (3)In this expression, the index t refers to the isoscalar ( t = 0 ) or isovector ( t = 1 ) channel.We refer to Refs. [22, 23, 24, 25, 26] for the definition of the densities ρ , τ , J , and J . Thereare 8 real-valued coupling constants C uu (cid:48) t for t = 0 , and uu (cid:48) = ( ρτ, J J, ρ ∆ ρ, ρ ∇ J ) . Thecase of C ρρt is a little different, since it is a function of the density ρ ( r ) , C ρρt = C ρρt + C ρρtD ρ γ ( r ) . (4)This density-dependent term is thus characterized by 5 parameters, the two C ρρt and C ρρt D and the exponent γ . Therefore, E nuc [ ρ ] is fully characterized by 13 parameters.The Coulomb term in (1) is computed at the Hartree-Fock approximation with theexchange term treated at the Slater approximation [24]. The pairing energy is computedfrom a surface-volume density-dependent pairing force V q ( r , r (cid:48) ) = V q (cid:20) − ρ ( r ) ρ c (cid:21) δ ( r − r (cid:48) ) , (5)where q here refers to the type of particle (proton or neutron) and ρ c = 0 . fm − .Including the pairing channel in the fit thus adds 2 more parameters to the fit. Ascustomary for zero-range pairing forces, a cut-off at E cut = 60 MeV is introduced tolimit the number of quasiparticles used when calculating the density matrix.In the context of fission, potential energy curves (or surfaces) are obtained byadding constraints on the density ρ . The constraints are typically the expectation valueof suitable operators on the hfb vacuum. In this work, we will consider only oneconstraint on the expectation value of the axial quadrupole moment q ≡ (cid:104) ˆ Q (cid:105) . Thetotal hfb energy at the deformation q is thus the scalar function E ( q ) . It implicitlydepends on the vector of coupling constants x ≡ { C uu (cid:48) t } of the energy density functional, E ≡ E ( q ; x ) . Following [27], we express the coupling constants C ρρt , C ρρtD , γ and C ρτt as a function of theparameters of infinite nuclear matter. Our analysis will be focused on the unedf1 HFB functional, for which the two tensor coupling constants are 0 and the vector effectivemass is set at the SLy4 value, m ∗ v = 1 . [28]. As a result, the actual vector x ofparameters is x = (cid:16) E NM , ρ sat , K NM , a NMsym , L
NMsym , m ∗ s , C ρ ∆ ρ , C ρ ∆ ρ , C ρ ∇ J , C ρ ∇ J , V n , V p (cid:17) (6)3LNL-JRNL-810569and the total number of parameters under consideration is n = 12 .Calculations were performed with the latest version of the code hfodd [29]. Thesolutions to the hfb equation are expanded on a deformed basis of N shell = 30 shells.For each value of the quadrupole moment, the oscillator length and axial quadrupoledeformation of the basis, which determines the ratios of oscillator frequencies ω z /ω x ,are set according to the empirical formula of [30]. The Gauss-Hermite integration meshcomprises 40 points in the x - and y -directions and 66 in the z -direction. We use linearconstraints on the quadrupole moment (and on the dipole moment to fix the positionof the center of mass), and the value of the Lagrange parameter is set at each iterationbased on the cranking approximation of the qrpa matrix.
3. Statistical methods
In this section we first introduce our Gaussian Process emulator of the potential energycurves E ( q ; x ) and then describe how this emulator will be used to generate a posteriordistribution of the parameters. To emulate the potential energy curves, we use a Gaussian Stochastic Process ( g a sp )emulator as described in [31]. To implement a g a sp , the data should be observed at thesame input locations for all curves, i.e., all potential energy curves must be defined onthe same grid of deformations q . A g a sp y ( · ) can be defined as y ( · ) ∼ g a sp ( µ ( · ) , σ C ( · , · )) , (7)with µ the mean function, σ the variance and C the correlation function. Here the meanfunction is modeled as a linear combination of basis functions h ( x ) ≡ ( h ( x ) , . . . , h q ( x )) and regression parameters κ ≡ ( κ , . . . , κ q ) so that µ ( x ) = E [ y ( x )] = h ( x ) · κ = q (cid:88) t =1 h t ( x ) κ t . (8)The correlation function is defined by the Matérn covariance structure. We apply acommon specific case referred to as the Matérn 5/2, C ( x , x (cid:48) ) = (cid:32) √ || x − x (cid:48) || γ + 5 || x − x (cid:48) || γ (cid:33) exp (cid:32) − √ || x − x (cid:48) || γ (cid:33) , (9)where γ is a hyper-parameter that is optimized to obtain a better fit. This correlationstructure is stationary, meaning that the correlation is assumed to only depend on thedistance between two observations, i.e., || x − x (cid:48) || , but not based on the location of thoseobservations. Note that each deformation q is treated independently of the others. Thismodel is fit using the RobustGaSP R package [32]. The default setting uses a constantmean for µ ( x ) but first subtracts the mean at each input location. The model is robustin terms of the estimation of hyperparameters such as γ . It avoids numerical issueswhile still achieving good predictive performance.4LNL-JRNL-810569 We now discuss the Bayesian framework that will be used to find posterior distributionsof the parameters. The philosophy of Bayesian inference starts with the fact that thereis a prior probability distribution representing the a priori beliefs or knowledge aboutthe parameter set. This prior state of knowledge is updated when presented with theobserved experimental data to give a posterior probability distribution representingthe current state of knowledge given the data. Formally, our goal is to determinethe posterior distribution of the parameters x represented as π ( x | y ) using observations y , . . . , y n using the formula π ( x | y ) ∝ L y ( x ) π ( x ) , (10)where π ( x ) is the prior distribution and L is the likelihood function of the data (thetheoretical prediction of the model) given the parameters. Here we use a use a normallikelihood function justified by assuming Gaussian measurement errors. Given a set { d i } i =1 ,n d of n d experimental values (in our case the three measurements in Table 2) forsome observables and { y i ( x ) } i =1 ,n d the theoretical predictions obtained from the g a sp emulator at x for these same observables, then the likelihood function reads L y ( x ) ∝ (cid:32) n d (cid:89) i =1 σ i (cid:33) nd exp (cid:32) − n d (cid:88) i =1 ( y i ( x ) − d i ) σ i (cid:33) , (11)where σ i is the variance for each experimental value. The advantage of Bayesianinference is that our result is a full distribution which allows us to make statements aboutthe uncertainty of our parameter estimates. The difficulty here is that an analytic form ofthe posterior distribution can only be found under specific distributional assumptions onthe prior and likelihood function. These assumptions are not always realistic, and as suchapproximations of the posterior are typically made using Markov Chain Monte Carlo( mcmc ) sampling techniques. These methods yield a large number of approximatelyindependent samples from a Markov chain that should have the same distribution asthe posterior. Here, we select non-informative prior distributions that are uniform overthe given parameter boundaries. This implies that we have no prior knowledge of whichregion of the parameter space is closer to the “real value”; we only assume that it willoccur somewhere in the specified region.We apply a Delayed Rejection Adaptive Metropolis mcmc algorithm [33]. TheDelayed Rejection algorithm samples a second value when a proposed value is rejected.This second proposal is typically drawn from a proposal distribution with a smallervariance. We use a maximum number of delayed rejections of one at each iteration.It is possible to do this multiple times, but it will greatly increase the computationtime of the mcmc algorithm if the maximum number of delayed rejections at eachiteration is too large. Delayed rejection is advantageous because it has been shown togive smaller asymptotic variances of the estimators from the chain [34]. The AdaptiveMetropolis component updates the covariance of the proposal distribution based onprevious iterations of the chain. This adapts the shape and size of the sampling5LNL-JRNL-810569 Table 1.
Range of variation around the unedf1
HFB value for the six couplingconstants included in this work. C ρ ∆ ρt and C ρ ∇ Jt are in MeV fm and V n and V p inMeV fm . C ρ ∆ ρ C ρ ∆ ρ V n V p C ρ ∇ J C ρ ∇ J ∆ x distribution and generally makes the mcmc more efficient. In our implementation thisupdate is only done during the burn-in phase.
4. Results
In this section, we focus on the one-dimensional potential energy curve of
Pu asas function of the axial quadrupole moment, from the ground state to the point ofscission. We first describe the calculations performed to generate the training data forthe Gaussian process analysis. We then discuss the emulator of the potential energycurve before determining the posterior distribution of the edf parameters.
Gaussian processes must be trained on some input data. In our case, these consist of aset of potential energy curves for a set of N parametrizations X = ( x , . . . , x N ) of theSkyrme functional. For each parametrization x k , we must calculate the full potentialenergy curve – here with triaxial deformations included. Since the parameter space of unedf1 HFB is 12-dimensional, the amount of calculation needed could quickly becomegigantic. We thus imposed two restrictions to alleviate the computational cost: • We reduced the number of varying coupling constants to 6 by keeping all bulkcoupling constants fixed at their unedf1
HFB value. This is justified since these arerelated to nuclear matter properties and the result of the optimization show thatthey are relatively well-constrained [35, 8]; • For each remaining coupling constant, we considered a rather small intervalof variation centered around the nominal unedf1
HFB value, I uu (cid:48) t = [ C uu (cid:48) t − ∆ C uu (cid:48) t , C uu (cid:48) t + ∆ C uu (cid:48) t ] . The quantities ∆ C uu (cid:48) t are listed in table 1 and correspondapproximately to the standard deviations of the unedf1 fit; see [35].The product set of all intervals I uu (cid:48) t defines an hypercube in the -dimensionalparameter space. We sampled this hypercube with a Latin Hypercube Sampling ( lhs )algorithm to determine the vector X of parameterizations for the training runs. Inpractice, as we will show below, N = 70 samples were sufficient to build a high-quality emulator of the model. These included 10 runs selected using the IntegratedMean Square Prediction Error criterion (IMSPE) [36] to improve estimation of thescission point. Figure 1 shows the deformation energy (with respect to the ground-state6LNL-JRNL-810569value) for each of the 70 samples for the case of Pu. The end-point of each curve,usually beyond q > b, represents the scission configuration: immediately beyondthis point, the total energy drops rapidly, which manifests itself by a discontinuity inthe energy curve. To increase the legibility of the figure, we did not represent thisdiscontinuity and simply stopped the curve at scission.
0 50 100 150 200 250 300 350 400Axial quadrupole moment q [b]-15-10 -5 0 5 10 D e f o r m a t i o n E n e r g y [ M e V ] Pu Figure 1.
Deformation energy curves in
Pu as a function of the axial quadrupolemoment q ≡ (cid:104) ˆ Q (cid:105) for a training set of 60 different parametrizations of the Skyrme edf . The rectangular grid guides the eye. In this section we fit the g a sp model to the potential energy curves. For the datasetof 70 training runs, the observations were not equally spaced, but they were dense.Thus, the points used to fit the model were selected by equal spacing via interpolation.Since the scission point was not needed for calibration, the emulator was only run upto q = 300 b. We are not only interested in having accurate estimation of the potentialenergy curve for the parameter sets we have observed, but also for any parameter setwithin the bounds of interest. Thus, we will use cross validation to look at the out ofsample error for our emulator, specifically applying leave-one-out cross validation. Thisis effective because it does not require additional samples, and shows what the modelpredictions would have been for each parameter set had the potential energy not been7LNL-JRNL-810569observed. Figure 2 gives the leave-one-out residuals for all 70 curves. To obtain leave-one-out residuals for curve y i , we first fit the model using the other 69 observations.Then, we use the predicted mean curve ˆ y i using the parameter set x i to obtain theresidual. −0.40.00.4 100 200 300 Quadropole Moment [b]
Lea v e − one − ou t r e s i dua l Figure 2.
Leave-one-out residuals (cid:15)
GaSP = E ( q ) − E GaSP ( q ) (in MeV) for the g a sp emulator as a function of the axial quadrupole moment. The results show that the largest errors (cid:15)
GaSP occur around the localmaxima/minima of the potential energy curve with the areas in between having muchsmaller residuals. In fact, a closer analysis of the emulation error shows that it is directlyconnected with the onset of triaxiality. This is better visualized in figure 3, which showsthe expectation values (cid:104) ˆ Q (cid:105) and (cid:104) ˆ Q (cid:105) as a function of the axial quadrupole moment.As a reminder, the degree of triaxiality of the nuclear shape is typically characterized bythe Bohr γ angle. With the conventions adopted for the multipole moments in hfodd ,we have γ = atan( (cid:104) ˆ Q (cid:105) / (cid:104) ˆ Q (cid:105) ) . Therefore, the shape is prolate axial if γ = 0 o (hence (cid:104) ˆ Q (cid:105) = 0 b), oblate axial if γ = 60 o ( (cid:104) ˆ Q (cid:105) = √ (cid:104) ˆ Q (cid:105) ), and it is maximally triaxial if γ = 30 o ( (cid:104) ˆ Q (cid:105) = (cid:112) / (cid:104) ˆ Q (cid:105) ).The onset of triaxiality at q ≈ b (depending on the parameter set), the returnto axial symmetry at q ≈ b, and the slightly trixial shapes between
110 b ≥ q ≥ b correspond exactly to the regions where the emulator error is maximal, reaching up to500 keV. This is most likely the consequence of spurious discontinuities that appear inthe self-consistent calculations when projecting on one-dimensional paths as in figure 2.Indeed, while the expectation value of the axial quadrupole moment q is constrained,8LNL-JRNL-810569every other moment of the nuclear surface is unconstrained. The variational principledictates that the hfb calculation will converge to the local minimum nearest to thestarting point. In the higher-dimensional space characterized by an additional collectivevariable q (cid:48) , one may find at point q = q two minima quasi-degenerate in energy butseparated by a barrier along the q (cid:48) direction: projecting such a surface on the q axiswould give a continuous, non-differentiable curve at q ; see discussion in [37]. This isexactly what seems to happen here, with q (cid:48) ≡ q controlling the degree of triaxiality.By contrast, the onset of octupole deformation around q ≈ b does not produce anynoticeable increase in emulation error.
20 40 60 80 100 120 140Quadrupole Moment Q [b] 0 2 4 6 8 10 12 14 M u l t i p o l e M o m e n t s ] Pu Plain: Q [b]Dashed: Q [b ] Figure 3.
Expectation value (cid:104) ˆ Q (cid:105) (triaxial quadrupole) and (cid:104) ˆ Q (cid:105) (axial octupole)of the multipole moments as a function of the axial quadrupole moment q ≡ (cid:104) ˆ Q (cid:105) . The size of this emulation error | (cid:15) GaSP | ≤ . MeV should be compared to otherestimates of relevant uncertainty: (i) for triaxial hfb calculations based on expansionsin the one-center harmonic oscillator basis, basis truncation errors alone can easily reacha few MeV [38] (ii) while the experimental value of the fission isomer excitation energyin
Pu is known to within ± . MeV [39], this number comes from a single experimentand older estimates differ by 0.6 MeV [40] (iii) the height of both fission barriers is amodel-dependent quantity extracted from fission cross sections, and the uncertainty istypically of the order of 1 MeV; see discussion in section 4.3 below.Figure 4 shows four examples of the centered potential energy curves along with9LNL-JRNL-810569the leave-one-out mean prediction and uncertainty from the g a sp emulator. Given the n = 70 samples, the centered potential energy curve E ( q, x i ) is obtained by subtractingat each point the mean value across all samples, E cent ( q, x i ) = E ( q, x i ) − n n (cid:88) k =1 E ( q, x k ) . (12)Centered potential energy curves with mostly negative values simply indicate that theoriginal potential energy at each deformation q is lower than the average value at thatdeformation over the n samples. The four parameter sets were chosen for visual clarityand are representative of the typical performance of the emulator. From this we can seethat the emulator is highly confident and also accurately captures the potential energycurves. There is more uncertainty at the minima and maxima of the curve, so the largerobserved deviations from the mean prediction still largely fall within the 95% credibleintervals. −4−202 100 200 300 Quadropole Moment [b] C en t e r ed D e f o r m a t i on E ne r g y [ M e V ] Parameter Set2123165
Figure 4.
Centered potential energy curves with g a sp predictions. Dashed linesrepresent the mean prediction of the g a sp model, and shaded areas represent 95%credible intervals. The solid lines represent the simulated values. Many relevant fission observables cannot be computed with a single collectivevariable q but require a higher-dimensional collective space. This is the case, for instance,with the distributions in charge, mass and kinetic energy of the fission fragments. Mostimportantly, these observables also need to be computed at scission, i.e., just before10LNL-JRNL-810569the system has split into two fragments. As mentioned in the previous section, in aone-dimensional space, scission corresponds to a point (the end-point of each curve infigure 1); in a collective space of dimensions D , scission configurations correspond to a ( D − -dimensional hypersurface. The characteristics of this surface are a function ofthe parameters x of the functional which we need to emulate. Before embarking in sucha project, which would imply running expensive calculations in a D -dimensional space,it is worth checking the ability of Gaussian processes to reproduce the scission pointalready in D = 1 . In practice, we try to emulate the location q scis of the discontinuity inthe potential energy curve as a function of x , as well as the number of particles Z H and N H of the heavy fragment at that point. These values are scalars and we use a simpleGaussian Process for prediction.The leave-one-out residuals are shown in figure 5 in the form of boxplots and “violin”plots. The upper and lower portions of the box portion of the boxplot show the 25thand 75th percentile of the data, and the horizontal line in the middle shows the median.The minimum and maximum values are indicated by the ends of the vertical line goingthrough the box; except for cases where values are far enough away form the box portionof the data (traditionally found as 1.5 × (75th percentile - 25th percentile) above orbelow the box portion). In these cases the individual points are plotted, such as thelargest negative residual for the scission point. The shaded area is the violin plot andit gives an estimate of the probability distribution of the data.Disregarding the outliers for the moment, we see that on average, the location ofthe scission point is reproduced within about 5 – 10 barns. The impact on the numberof particles in the heavy fragment is of the order of half a particle in total. From aphysics perspective, these numbers may look small, but one should bear in mind thatthey were obtained for the case of a one-dimensional collective space only: they arelikely to increase with the size of the collective space and/or the size of the parameterspace that we try to cover with the emulator even as they decrease with the numberof training runs. It is not unreasonable to believe that such emulators could lead touncertainties σ A > in the number of particles at scission. Such uncertainties have tobe added to the others related to the very concept of scission; see discussion in [3]. It isalso interesting to note that even though there is only one outlier for the emulation of thescission point, there are 3 of them for Z H and A H . This could be another manifestationof the aforementioned discontinuities, this time at scission. This effect is also likelyto be further magnified when moving to higher-dimensional collective spaces. One mayconclude from these results that Gaussian processes may not be the best tool to emulatethe characteristics of scission. With the relative accuracy of the emulator established, we now proceed to obtainposterior distributions of the parameters using the g a sp emulator by conditioning on theexperimental value of fission isomer excitation energy and fission barriers. This exercise11LNL-JRNL-810569 −20−10010 Scission Point
Lea v e − O ne − O u t R e s i dua l s −1.00−0.75−0.50−0.250.000.25 ZH −2−10 AH Figure 5.
Boxplots overlaying violin plots for the leave-one-out prediction error forthe location of the scission point, the number of protons in the heavy fragment (ZH)and the number of particles in the heavy fragment (AH)
Table 2.
Summary of experimental and empirical information about the fissionisomer excitation energy of
Pu and its two fission barriers. The last column is ourweighted average. All units are MeV. Data for Bjornholm & Lynn is from [41]; forCapote from [42]; for Hilaire from [43]; for Hunyadi from [39] and for Singh from [40].Bjornholm Capote Hilaire Hunyadi Singh Weighted E FI ± ± ± E A ± E B ± presents two difficulties. First, there is little experimental information on the excitationenergy of fission isomers and, more generally, the band-head of the lowest rotationalband built on superdeformed minima [40]. Second, fission barriers are extracted fromfission cross-sections in a model-dependent procedure. Typically, this analysis is basedon assuming a one-dimensional, inverted parabola for the barrier [41]. As a result,fission barriers are not genuine experimental data, but empirical one, sometimes called“metadata”. Table 2 summarizes what is available in the literature for the case of Pu.To generate the posterior distribution we first need to define our likelihood function.When using the normal likelihood (11), posterior distributions generated about 10% of12LNL-JRNL-810569
Table 3.
Values of the the six coupling constants included in this work for theoriginal unedf1
HFB parameterization and the map estimate. C ρ ∆ ρt and C ρ ∇ Jt are inMeV fm and V n and V p in MeV fm . C ρ ∆ ρ C ρ ∆ ρ V n V p C ρ ∇ J C ρ ∇ J unedf1 HFB -45.600 -143.935 -234.380 -260.437 -73.946 -51.913 map -46.157 -139.972 -245.287 -250.964 -77.353 -70.143 potential energy curves with E A < E B . Since there is a relative consensus in the physicscommunity that E A > E B , we modified the likelihood in (11) to be a truncated normallikelihood: L ( y | x ) is set to be extremely small if E A < E B so that the corresponding x will not be accepted during mcmc iterations.The values of d i for the fission isomer excitation energy and the two fission barrierswere calculated as a weighted average of the experimental or empirical values listedin table 2 based on subjective confidence estimates. Specifically, for E FI , the Hunyadivalue accounted for 60% of the “experimental” value while the Bjornholm and Singhnumbers accounted for 30% and 10%, respectively. For both fission barriers, theCapote value accounted for 50% and the Hilaire and Bjornholm numbers for 25% each.The last column of table 2 gives the actual weighted values d i that we used in theregression. For the calculation of the likelihood (11), the variances σ i were taken as σ i = (cid:113) σ + σ i, g a sp , where σ exp is the experimental error, see table 4, and σ i, g a sp , isthe sum of the emulator standard deviation at the two locations used to calculate eachexcitation energy energy.To generate the posterior samples, the first 10,000 mcmc samples were discardedin the burn-in period. The chain ran for a total of 7.5 million samples, but was thinneddown to 100,000 to ensure the posterior samples were uncorrelated. The maximum aposteriori ( map ) estimator x ∗ was found by taking the parameter set with the highestposterior probability. Table 3 compares the values of the six coupling constants forthe map estimate and for the original unedf1 HFB parametrization. Even though theintervals of variation around the unedf1
HFB values was relatively small as shown intable 1, most parameters have changed significantly. The two-dimensional bivariaterepresentation of the full 6-dimensional posterior distribution is shown in figure 6.Clearly, the local fit pins down the value of C ρ ∆ ρ and V n , which is compatible withthe analysis of [44, 45] suggesting that surface properties, which are largely dependenton the interplay between C ρ ∆ ρ and V n , are highly correlated with the minima andmaxima of potential energy curves. Conversely, all other coupling constants end attheir boundaries, which indicates that they are not well-constrained by deformationproperties.Given the x ∗ parameterization for the map estimate, we can then run the g a sp emulator for this set of coupling constants and reconstruct the entire potential energycurve. The results are shown in figure 7 together with the curve obtained from the13LNL-JRNL-810569 Figure 6.
Posterior estimates of the parameters. The blue lines on the diagonal andthe white dots with black outlines on the off diagonals indicate the map values. initial unedf1
HFB parameterization. The shaded band around the map curve showsthe 95% uncertainty from the g a sp emulator. Even if we only calibrated 3 excitationenergies, E FI , E A and E B , the Bayesian regression modifies the entire deformation energycurve. This could have a significant effect on observables such as spontaneous fission halflives τ SF which are related to the quantum-mechanical tunneling probability through thebarrier and thus depend both on the barrier height but also its width. Qualitatively, the τ SF for Pu computed from the map estimate could thus be a few orders of magnitudesmaller than the value obtained with the initial parameterization of the functional, sincethe barriers are both lower and narrower.The excitation energy and fission barrier heights can be extracted easily from thepotential energy curves and are listed in table 4. The uncertainties are from the emulatoruncertainty at the ground state and the relevant location for each energy. We recall thatthe empirical value used to determine the map solution is a weighted average of threedifferent numbers; see table 2. In the case of the second barrier, one of the numbers issignificantly different from the others, which leads to a larger uncertainty and explainswhy the map result is further away from the data.14LNL-JRNL-810569 −2.50.02.55.0 100 200 300
Quadropole Moment [b] D e f o r m a t i on E ne r g y [ M e V ] Parameter SetMAPUNEDF1
HFB
MAP Estimate
Figure 7.
Potential energy curve generated using the map and unedf1
HFB parameterizations. In both cases, the curve was obtained from the g a sp emulator,which provides the 95% confidence region. Table 4.
Energy of the fission isomer ( E FI ), height of the first ( E A ) and second( E B ) barrier in MeV for the original unedf1 HFB functional, the parameterization thatmaximizes the likelihood and empirical data. We also indicate the experimental error( σ exp ) that was used in the regression. unedf1 HFB map σ map Empirical σ exp E FI . E A . E B .
5. Conclusions
In this work, we investigated the use of Gaussian processes to quantify uncertainties innuclear deformation properties. We built an emulator of the entire one-dimensionalpotential energy curve in
Pu that is valid for a small 6-dimensional hypercubearound the unedf1
HFB parameterization of the Skyrme functional. When the potentialenergy curve follows a fully-connected path in the (infinite-dimensional) collectivespace, the numerical precision of the emulator is excellent (less than 100 keV); when“discontinuities” occur as, e.g., triggered by the onset of triaxiality, the quality of theemulation is a little degraded but remains rather good (less than 500 keV). The location15LNL-JRNL-810569of the scission point, which is marked as a discontinuity in the potential energy, is harderto pin-down, with uncertainties of the order of ∆ q ≈ ± b. We used our emulator todetermine the posterior distribution of the coupling constants of the Skyrme energyfunctional conditioned on the excitation energy of the fission isomer and the height ofthe two fission barriers. Most of the uncertainty comes from the empirical data.Gaussian processes belong to the class of supervised learning techniques, whichthemselves are part of the broader field of machine learning. In nuclear density functionaltheory, these techniques can be applied on at least three different classes of “data”: (1)the observable of interest, i.e., the total energy, radius, cross-section, etc. (2) matrixelements, e.g., of the mean field, pairing field or of some collective Hamiltonian and (3)the single-particle or quasiparticle wave functions. This work showed that there maybe somewhat limited use in trying to emulate directly the observable of interest. Therecent work of [46] suggests that working with matrix elements may hold more promise.Ultimately, one may have to directly work with individual wave functions themselves. Acknowledgment
This work was partly performed under the auspices of the U.S. Department of Energyby Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344and was supported by the LLNL-LDRD Program under Project No. 18-ERD-008.Computing support also came from the Lawrence Livermore National Laboratory(LLNL) Institutional Computing Grand Challenge program.
References [1] M. R. Mumpower, R. Surman, G. C. McLaughlin, and A. Aprahamian. The impact of individualnuclear properties on r-process nucleosynthesis.
Prog. Part. Nucl. Phys. , 86:86, 2016.[2] S. A. Giuliani, Z. Matheson, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sadhukhan,B. Schuetrumpf, N. Schunck, and P. Schwerdtfeger. Colloquium: Superheavy elements:Oganesson and beyond.
Rev. Mod. Phys. , 91(1):011001, 2019.[3] N Schunck and L M Robledo. Microscopic theory of nuclear fission: A review.
Rep. Prog. Phys. ,79(11):116301, 2016.[4] R.M. Dreizler and E.K.U Gross.
Density Functional Theory: An Approach to the Quantum Many-Body Problem . Springer-Verlag, 1990.[5] R. Eschrig.
Fundamentals of Density Functional Theory . Teubner, Leipzig, 1996.[6] Nicolas Schunck.
Energy Density Functional Methods for Atomic Nuclei.
IOP Expanding Physics.IOP Publishing, Bristol, UK, 2019. OCLC: 1034572493.[7] N. Schunck, J. D. McDonnell, D. Higdon, J. Sarich, and S. Wild. Quantification of Uncertaintiesin Nuclear Density Functional Theory.
Nucl. Data Sheets , 123:115, 2015.[8] Nicolas Schunck, Jordan D. McDonnell, Jason Sarich, Stefan M. Wild, and Dave Higdon. Erroranalysis in nuclear density functional theory.
J. Phys. G: Nucl. Part. Phys. , 42(3):034024, 2015.[9] N. Schunck, J. D. McDonnell, D. Higdon, J. Sarich, and S. M. Wild. Uncertainty Quantificationand Propagation in Nuclear Density Functional Theory.
Eur. Phys. J. A , 51(12):1, 2015.[10] Y. Gao, J. Dobaczewski, M. Kortelainen, J. Toivanen, and D. Tarpanov. Propagation ofuncertainties in the Skyrme energy-density-functional model.
Phys. Rev. C , 87(3):034324, 2013. [11] S. Goriely and R. Capote. Uncertainties of mass extrapolations in Hartree-Fock-Bogoliubov massmodels.
Phys. Rev. C , 89(5):054318, 2014.[12] R. Utama, J. Piekarewicz, and H. B. Prosper. Nuclear mass predictions for the crustal compositionof neutron stars: A Bayesian neural network approach.
Phys. Rev. C , 93(1):014311, 2016.[13] T. Haverinen and M. Kortelainen. Uncertainty propagation within the UNEDF models.
J. Phys.G: Nucl. Part. Phys. , 44(4):044008, 2017.[14] Léo Neufcourt, Yuchen Cao, Witold Nazarewicz, and Frederi Viens. Bayesian approach to model-based extrapolation of nuclear observables.
Phys. Rev. C , 98(3):034318, 2018.[15] Léo Neufcourt, Yuchen Cao, Witold Nazarewicz, Erik Olsen, and Frederi Viens. Neutron DripLine in the Ca Region from Bayesian Model Averaging.
Phys. Rev. Lett. , 122(6):062502, 2019.[16] M. Kortelainen, J. Erler, W. Nazarewicz, N. Birge, Y. Gao, and E. Olsen. Neutron-skinuncertainties of Skyrme energy density functionals.
Phys. Rev. C , 88(3):031305, 2013.[17] P.-G. Reinhard and W. Nazarewicz. Information content of the low-energy electric dipole strength:Correlation analysis.
Phys. Rev. C , 87(1):014324, 2013.[18] N. Paar, Ch. C. Moustakidis, T. Marketin, D. Vretenar, and G. A. Lalazissis. Neutron starstructure and collective excitations of finite nuclei.
Phys. Rev. C , 90(1):011304, 2014.[19] J. D. McDonnell, N. Schunck, D. Higdon, J. Sarich, S. M. Wild, and W. Nazarewicz. UncertaintyQuantification for Nuclear Density Functional Theory and Information Content of NewMeasurements.
Phys. Rev. Lett. , 114(12):122501, 2015.[20] J.-P. Blaizot and G. Ripka.
Quantum Theory of Finite Systems . The MIT Press, Cambridge,1985.[21] P. Ring and P. Schuck.
The Nuclear Many-Body Problem . Texts and Monographs in Physics.Springer, 2004.[22] Y. M. Engel, D. M. Brink, K. Goeke, S. J. Krieger, and D. Vautherin. Time-dependent Hartree-Fock theory with Skyrme’s interaction.
Nucl. Phys. A , 249(2):215, 1975.[23] J. Dobaczewski and J. Dudek. Time-Odd Components in the Rotating Mean Field and IdenticalBands.
Acta Phys. Pol. B , 27(1):45, 1996.[24] Michael Bender, Paul-Henri Heenen, and Paul-Gerhard Reinhard. Self-consistent mean-fieldmodels for nuclear structure.
Rev. Mod. Phys. , 75(1):121, 2003.[25] E. Perlińska, S. G. Rohoziński, J. Dobaczewski, and W. Nazarewicz. Local density approximationfor proton-neutron pairing correlations: Formalism.
Phys. Rev. C , 69(1):014316, 2004.[26] T. Lesinski, M. Bender, K. Bennaceur, T. Duguet, and J. Meyer. Tensor part of the Skyrmeenergy density functional: Spherical nuclei.
Phys. Rev. C , 76(1):014312, 2007.[27] M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, andS. Wild. Nuclear energy density optimization.
Phys. Rev. C , 82(2):024313, 2010.[28] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer. A Skyrme parametrizationfrom subnuclear to neutron star densities Part II. Nuclei far from stabilities.
Nucl. Phys. A ,635(1):231, 1998.[29] N. Schunck, J. Dobaczewski, W. Satuła, P. Bączyk, J. Dudek, Y. Gao, M. Konieczka, K. Sato,Y. Shi, X. B. Wang, and T. R. Werner. Solution of the Skyrme-Hartree–Fock–Bogolyubovequations in the Cartesian deformed harmonic-oscillator basis. (VIII) HFODD (v2.73y): A newversion of the program.
Comput. Phys. Commun. , 216:145, 2017.[30] N. Schunck, D. Duke, H. Carr, and A. Knoll. Description of induced nuclear fission with Skyrmeenergy functionals: Static potential energy surfaces and fission fragment properties.
Phys. Rev.C , 90(5):054305, 2014.[31] Mengyang Gu, Xiaojing Wang, James O Berger, et al. Robust gaussian stochastic processemulation.
The Annals of Statistics , 46(6A):3038–3066, 2018.[32] Mengyang Gu, Jesus Palomo, and James Berger.
RobustGaSP: Robust Gaussian Stochastic ProcessEmulation , 2019. R package version 0.5.7.[33] Heikki Haario, Marko Laine, Antonietta Mira, and Eero Saksman. Dram: efficient adaptive mcmc.
Statistics and computing , 16(4):339–354, 2006. [34] Luke Tierney and Antonietta Mira. Some adaptive monte carlo methods for bayesian inference.
Statistics in medicine , 18(17-18):2507–2515, 1999.[35] M. Kortelainen, J. McDonnell, W. Nazarewicz, P.-G. Reinhard, J. Sarich, N. Schunck, M. V.Stoitsov, and S. M. Wild. Nuclear energy density optimization: Large deformations.
Phys.Rev. C , 85(2):024304, 2012.[36] Thomas J Santner, Brian J Williams, William Notz, and Brain J Williams.
The design andanalysis of computer experiments , volume 1. Springer, 2003.[37] N. Dubray and D. Regnier. Numerical search of discontinuities in self-consistent potential energysurfaces.
Comput. Phys. Commun. , 183(10):2035, 2012.[38] N. Schunck. Density Functional Theory Approach to Nuclear Fission.
Acta Phys. Pol. B , 44:263,2013.[39] M. Hunyadi, D. Gassmann, A. Krasznahorkay, D. Habs, P. G. Thirolf, M. Csatlós, Y. Eisermann,T. Faestermann, G. Graw, J. Gulyás, R. Hertenberger, H. J. Maier, Z. Máté, A. Metz, and M. J.Chromik. Excited superdeformed K π =0+ rotational bands in β -vibrational fission resonancesof Pu.
Phys. Lett. B , 505(1):27–35, 2001.[40] Balraj Singh, Roy Zywina, and Richard B. Firestone. Table of superdeformed nuclear bands andfission isomers: (october 2002).
Nucl. Data Sheets , 97(2):241, 2002.[41] S. Bjørnholm and J. E. Lynn. The double-humped fission barrier.
Rev. Mod. Phys. , 52(4):725,1980.[42] R. Capote, M. Herman, P. Obložinský, P.G. Young, S. Goriely, T. Belgya, A.V. Ignatyuk, A.J.Koning, S. Hilaire, V.A. Plujko, M. Avrigeanu, O. Bersillon, M.B. Chadwick, T. Fukahori,Zhigang Ge, Yinlu Han, S. Kailas, J. Kopecky, V.M. Maslov, G. Reffo, M. Sin, E.Sh.Soukhovitskii, and P. Talou. RIPL – Reference Input Parameter Library for Calculation ofNuclear Reactions and Nuclear Data Evaluations.
Nucl. Data Sheets , 110(12):3107, 2009.[43] S. Goriely, S. Hilaire, A. J. Koning, M. Sin, and R. Capote. Towards a prediction of fission crosssections on the basis of microscopic nuclear inputs.
Phys. Rev. C , 79(2):024612, 2009.[44] N. Nikolov, N. Schunck, W. Nazarewicz, M. Bender, and J. Pei. Surface symmetry energy ofnuclear energy density functionals.
Phys. Rev. C , 83(3):034305, 2011.[45] W. Ryssens, M. Bender, K. Bennaceur, P.-H. Heenen, and J. Meyer. Impact of the surface energycoefficient on the deformation properties of atomic nuclei as predicted by Skyrme energy densityfunctionals.
Phys. Rev. C , 99(4):044315, 2019.[46] Raphaël-David Lasseri, David Regnier, Jean-Paul Ebran, and Antonin Penon. Taming NuclearComplexity with a Committee of Multilayer Neural Networks.
Phys. Rev. Lett. , 124(16):162502,2020., 124(16):162502,2020.