Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

A universal reduced basis for the calibration of covariant energy density functionals

Amy L. Anderson aanderson6@fsu.edu    J. Piekarewicz jpiekarewicz@fsu.edu Department of Physics, Florida State University, Tallahassee, FL 32306, USA
(June 3, 2024)
Abstract

The reduced basis method is used to construct a “universal” basis of Dirac orbitals that may be applicable throughout the nuclear chart to calibrate covariant energy density functionals. Relative to our earlier work using the non-relativistic Schrödinger equation, the Dirac equation adds an extra layer of complexity due to the existence of negative energy states. However, once this problem is mitigated, the resulting reduced basis is able to accurately and efficiently reproduce the high-fidelity model at a fraction of the computational cost. We are confident that the resulting reduced basis will serve as a foundational element in developing rapid and accurate emulators. In turn, these emulators will play a critical role in the Bayesian optimization of covariant energy density functionals.

I Introduction

Since first suggested by the editors of the Physical Review in 2011 PRA-Editors (2011), the theoretical nuclear physics community has embraced the notion of including uncertainty estimates in papers involving theoretical calculations of physical observables. To properly quantify theoretical uncertainties and explore correlations among physical observables, some of the initial efforts were devoted to the accurate calibration of a variety of nuclear models Reinhard and Nazarewicz (2010); Kortelainen et al. (2010); Fattoyev and Piekarewicz (2011); Reinhard et al. (2013); Dobaczewski et al. (2014); Piekarewicz et al. (2015); Chen and Piekarewicz (2014); Wesolowski et al. (2016); Yuksel et al. (2019); Drischler et al. (2020a, b); Furnstahl et al. (2020). Further, some earlier applications of artificial neural networks used to predict nuclear masses and charge radii that lie beyond the experimental reach Athanassopoulos et al. (2004); Akkoyun et al. (2013); Bayram et al. (2014); Gernoth et al. (1993), have now been extended to include Bayesian statistics for the proper quantification of theoretical uncertainties Utama et al. (2016a, b); Utama and Piekarewicz (2018); Neufcourt et al. (2019, 2020); Lovell et al. (2022); Saito et al. (2024). Finally, within the last few years, the Reduced Basis Method (RMB) have been incorporated into the arsenal of machine-learning tools dedicated to accelerate calculations of complex nuclear systems Frame et al. (2018); König et al. (2020); Drischler et al. (2021); Bonilla et al. (2022); Melendez et al. (2022); Giuliani et al. (2023); Anderson et al. (2022); Odell et al. (2024); Godbey et al. .

Perhaps the main virtue of the RBM is the ability to create a powerful framework for building efficient and accurate emulators. Even with advances in computer power, algorithmic developments, and physical insights, a full Bayesian calibration of nuclear models is generally unattainable without the help of emulators. Besides providing robust statistical uncertainties, Bayesian inference is often used to identify strong correlations between a property that is inaccessible and a surrogate observable that may be determined experimentally. However, Bayesian inferences that rely on Monte Carlo methods require a large number of evaluations of the same set of observables for many different realizations of the model parameters. This can be computationally expensive and highly impractical. Instead, emulators transform such a computationally-intractable problem into a low-dimensional calculation, thereby allowing enormous gains in computational speed with a minimal loss in precision. In essence, the reduced basis method, which falls under the general rubric of “Reduced Order Models” Brunton and Kutz (2019), encapsulate a set of dimensionality reduction approaches that aim at speeding up computations by approximating the solution of the high-fidelity problem by retaining only a handful of the most important components, namely, the coefficients of the reduced basis. The core assumption behind the RBM is that the solutions of the problem of interest vary smoothly across a subspace defined by the model parameters. Thus, the solution for an arbitrary set of parameters can be accurately approximated by a linear combination of exact solutions obtained for a small set of parameters.

Unlike other approaches, the RBM does not rely on uncontrolled extrapolations that aim to predict the behavior of the system outside its range of validity. Rather, the reduced basis provides a hierarchical scheme to select the most important basis states. As such, the most useful part of the RBM is the “R”, that instead of using the entire Hilbert space, selects an optimal basis for describing the eigenstates of a Hamiltonian with an arbitrary set of model parameters with a much smaller (reduced) set of functions. While the RBM have only recently been embraced by the nuclear physics community Bonilla et al. (2022); Melendez et al. (2022), it is a mature field that facilitates the solution of challenging numerical problems in a range of different areas Benner et al. (2020).

In this paper we will use eigenvector continuation Frame et al. (2018), a particular implementation of the RBM Quarteroni et al. (2015); Hesthaven et al. (2016), that has been shown to be highly successful for calculating the eigenvalues and eigenvectors of a Hamiltonian matrix—even when the “training” set lies far away from the physically relevant model parameters Frame et al. (2018). While a direct matrix diagonalization in the entire Hilbert space is computationally expensive, the RBM can generate the same low-energy states with minimal computational burden. Several areas in nuclear physics have benefited from the development of efficient and accurate emulators using the RBM König et al. (2020); Furnstahl et al. (2020); Drischler et al. (2021). Indeed, in a recent work we have constructed a “universal” reduced basis for the solution of a scaled Schrödinger equation that successfully reproduced the single-particle spectrum of a variety of nuclei Anderson et al. (2022). Ultimately, however, we aim to create robust emulators for the calibration of covariant energy density functionals Giuliani et al. (2023). Hence, the extension to the relativistic domain now demands the solution of the Dirac equation. Creating a “universal” basis of Dirac orbitals to compute the single-particle structure of a variety of nuclei is the main goal of this work.

The main justification behind this goal is twofold. First, we use the independent particle model, one of the fundamental pillars of nuclear structure. As argued by Bohr and Mottelson: “the relatively long mean free path of the nucleons implies that the interactions primarily contribute a smoothly varying average potential in which the particles move independently” Bohr and Mottelson (1998). While other important effects, such as pairing correlations, require a modification to the independent particle paradigm, we avoid these complications by focusing exclusively on the bulk properties of magic and semi-magic nuclei. Indeed, this kind of nuclei have been the bedrock of our calibration protocol; see for example Todd-Rutel and Piekarewicz (2005); Fattoyev et al. (2010); Chen and Piekarewicz (2014, 2015). Second, Density Functional Theory (DFT) guarantees that functional minimization of a “suitable” energy density functional yields both the exact ground-state energy and one-body density of the complicated many-body system Hohenberg and Kohn (1964). However, since DFT offers no guidance on how to construct such a suitable energy density functional, Kohn and Sham replaced the complex interacting system by an equivalent system of non-interacting fermions moving in an average potential rich enough to capture all the essential physics Kohn and Sham (1965). In the next sections we demonstrate how a physically-motivated average potential may be used to generate a spherical basis of Dirac orbitals that could be used through the entire nuclear chart.

The paper has been organized as follows. In Sec.II we review the structure of the Dirac equation and develop the formalism necessary to generate a reduced basis of Dirac orbitals. We then test in Sec.III the performance of the RBM in reproducing the high-fidelity results. Finally, we conclude in Sec.IV with a short summary and a vision for future work.

II Formalism

In 1928, Paul Dirac proposed an alternative to the Schrödinger equation for the description of electrons by incorporating the principles of special relativity. One of the most remarkable consequences of demanding that time and space be treated on equal footing is the unavoidable emergence of negative energy states Dirac (1982). As we will show below, the appearance of negative energy states requires special attention when constructing a robust reduced basis.

II.1 Dirac equation for spherically symmetric potentials

As mentioned above, we aim to build a RBM that serves as input for Bayesian optimization of covariant energy density functionals informed by the ground state properties of magic and semi-magic nuclei. As such, the training points for building the reduced basis only require the solution of the Dirac equation for spherically symmetric potentials with a Hamiltonian of the following form:

H=𝜶𝐏+β(MS(r))+V(r),𝐻𝜶𝐏𝛽𝑀𝑆𝑟𝑉𝑟H={\bm{\alpha}}\cdot{\bf P}+\beta\Big{(}M-S(r)\Big{)}+V(r),italic_H = bold_italic_α ⋅ bold_P + italic_β ( italic_M - italic_S ( italic_r ) ) + italic_V ( italic_r ) , (1)

where 𝜶𝜶{\bm{\alpha}}bold_italic_α and β𝛽\betaitalic_β are Dirac matrices, 𝐏=i𝐏𝑖bold-∇{\bf P}\!=\!-i{\bm{\nabla}}bold_P = - italic_i bold_∇ is the momentum operator, and M𝑀Mitalic_M is the bare nucleon mass. In turn, S(r)𝑆𝑟S(r)italic_S ( italic_r ) and V(r)𝑉𝑟V(r)italic_V ( italic_r ) are spherically symmetric Lorentz scalar and time-like vector potentials that underpin the successful relativistic mean field theory (RMF). The hallmark of the RMF framework is strong scalar and time-like vector potentials that cancel to produce a relatively weak central potential but which add coherently to account for the strong spin-orbit splitting observed in atomic nuclei Horowitz and Serot (1981); Serot and Walecka (1986).

The spherical symmetry of the Hamiltonian guarantees that the eigenstates of the Dirac Hamiltonian can be classified according to their total angular momentum j𝑗jitalic_j and its projection m𝑚mitalic_mSakurai (1967). That is,

𝒰nκm(𝐫)=1r(gnκ(r)𝒴+κm(𝐫^)ifnκ(r)𝒴κm(𝐫^)),subscript𝒰𝑛𝜅𝑚𝐫1𝑟subscript𝑔𝑛𝜅𝑟subscript𝒴𝜅𝑚^𝐫𝑖subscript𝑓𝑛𝜅𝑟subscript𝒴𝜅𝑚^𝐫{\cal U}_{n\kappa m}({\bf r})=\frac{1}{r}\left(\begin{array}[]{c}\phantom{i}g_% {n\kappa}(r){\cal Y}_{+\kappa\mbox{}m}(\hat{\bf r})\\ if_{n\kappa}(r){\cal Y}_{-\kappa\mbox{}m}(\hat{\bf r})\end{array}\right),caligraphic_U start_POSTSUBSCRIPT italic_n italic_κ italic_m end_POSTSUBSCRIPT ( bold_r ) = divide start_ARG 1 end_ARG start_ARG italic_r end_ARG ( start_ARRAY start_ROW start_CELL italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_r ) caligraphic_Y start_POSTSUBSCRIPT + italic_κ italic_m end_POSTSUBSCRIPT ( over^ start_ARG bold_r end_ARG ) end_CELL end_ROW start_ROW start_CELL italic_i italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_r ) caligraphic_Y start_POSTSUBSCRIPT - italic_κ italic_m end_POSTSUBSCRIPT ( over^ start_ARG bold_r end_ARG ) end_CELL end_ROW end_ARRAY ) , (2)

where n𝑛nitalic_n is the principal quantum number and the spin spherical harmonics 𝒴κmsubscript𝒴𝜅𝑚{\cal Y}_{\kappa\mbox{}m}caligraphic_Y start_POSTSUBSCRIPT italic_κ italic_m end_POSTSUBSCRIPT are obtained by coupling the orbital angular l𝑙litalic_l momentum to the intrinsic nucleon spin-1/2 to obtain a total angular momentum j𝑗jitalic_j. That is,

|κm=|l12jm.ket𝜅𝑚ket𝑙12𝑗𝑚|\kappa\,m\rangle=\Big{|}l\frac{1}{2}jm\Big{\rangle}.| italic_κ italic_m ⟩ = | italic_l divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_j italic_m ⟩ . (3)

In particular, the orbital angular momentum l𝑙litalic_l and the total angular momentum j𝑗jitalic_j are obtained from the generalized angular momentum κ𝜅\kappaitalic_κ as follows:

j=|κ|12;l={κ,ifκ>0;(1+κ),ifκ<0.formulae-sequence𝑗𝜅12𝑙cases𝜅if𝜅01𝜅if𝜅0j=|\kappa|\!-\!\frac{1}{2}\;;\quad l=\begin{cases}\kappa\;,&{\rm if}\;\kappa>0% ;\\ -(1+\kappa)\;,&{\rm if}\;\kappa<0.\end{cases}italic_j = | italic_κ | - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ; italic_l = { start_ROW start_CELL italic_κ , end_CELL start_CELL roman_if italic_κ > 0 ; end_CELL end_ROW start_ROW start_CELL - ( 1 + italic_κ ) , end_CELL start_CELL roman_if italic_κ < 0 . end_CELL end_ROW (4)

Unlike the case of the spherically symmetric Schrödinger equation, the orbital angular momentum is not a good quantum number. Instead, the orbital angular momentum of the upper and lower components differ by one unit.

By exploiting the spherical symmetry of the problem, one can eliminate the entire angular dependence and reduce the problem to a set of first order, coupled differential equations:

(ddr+κr)gnκ(r)[E+MS(r)V(r)]fnκ(r)=0,𝑑𝑑𝑟𝜅𝑟subscript𝑔𝑛𝜅𝑟delimited-[]𝐸𝑀𝑆𝑟𝑉𝑟subscript𝑓𝑛𝜅𝑟0\displaystyle\left(\frac{d}{dr}+\frac{\kappa}{r}\right)g_{n\kappa}(r)-\left[E+% M-S(r)-V(r)\right]f_{n\kappa}(r)=0,( divide start_ARG italic_d end_ARG start_ARG italic_d italic_r end_ARG + divide start_ARG italic_κ end_ARG start_ARG italic_r end_ARG ) italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_r ) - [ italic_E + italic_M - italic_S ( italic_r ) - italic_V ( italic_r ) ] italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_r ) = 0 , (5a)
(ddrκr)fnκ(r)+[EM+S(r)V(r)]gnκ(r)=0.𝑑𝑑𝑟𝜅𝑟subscript𝑓𝑛𝜅𝑟delimited-[]𝐸𝑀𝑆𝑟𝑉𝑟subscript𝑔𝑛𝜅𝑟0\displaystyle\left(\frac{d}{dr}-\frac{\kappa}{r}\right)f_{n\kappa}(r)+\left[E-% M+S(r)-V(r)\right]g_{n\kappa}(r)=0.( divide start_ARG italic_d end_ARG start_ARG italic_d italic_r end_ARG - divide start_ARG italic_κ end_ARG start_ARG italic_r end_ARG ) italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_r ) + [ italic_E - italic_M + italic_S ( italic_r ) - italic_V ( italic_r ) ] italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_r ) = 0 . (5b)

In search of a universal basis, it is convenient to rescale the above equations by an intrinsic distance scale cA1/3proportional-to𝑐superscript𝐴13c\!\propto A^{1/3}italic_c ∝ italic_A start_POSTSUPERSCRIPT 1 / 3 end_POSTSUPERSCRIPT characteristic of each individual nucleus. Hence, by defining the dimensionless distance as xr/c𝑥𝑟𝑐x\!\equiv\!r/citalic_x ≡ italic_r / italic_c, the coupled Dirac equation may now be written as

(ddx+κx)gnκ(x)[ϵ+μσ(x)ω(x)]fnκ(x)=0,𝑑𝑑𝑥𝜅𝑥subscript𝑔𝑛𝜅𝑥delimited-[]italic-ϵ𝜇𝜎𝑥𝜔𝑥subscript𝑓𝑛𝜅𝑥0\displaystyle\left(\frac{d}{dx}+\frac{\kappa}{x}\right)g_{n\kappa}(x)-\left[% \epsilon+\mu-\sigma(x)-\omega(x)\right]f_{n\kappa}(x)=0,( divide start_ARG italic_d end_ARG start_ARG italic_d italic_x end_ARG + divide start_ARG italic_κ end_ARG start_ARG italic_x end_ARG ) italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) - [ italic_ϵ + italic_μ - italic_σ ( italic_x ) - italic_ω ( italic_x ) ] italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) = 0 , (6a)
(ddxκx)fnκ(x)+[ϵμ+σ(x)ω(x)]gnκ(x)=0,𝑑𝑑𝑥𝜅𝑥subscript𝑓𝑛𝜅𝑥delimited-[]italic-ϵ𝜇𝜎𝑥𝜔𝑥subscript𝑔𝑛𝜅𝑥0\displaystyle\left(\frac{d}{dx}-\frac{\kappa}{x}\right)f_{n\kappa}(x)+\left[% \epsilon-\mu+\sigma(x)-\omega(x)\right]g_{n\kappa}(x)=0,( divide start_ARG italic_d end_ARG start_ARG italic_d italic_x end_ARG - divide start_ARG italic_κ end_ARG start_ARG italic_x end_ARG ) italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) + [ italic_ϵ - italic_μ + italic_σ ( italic_x ) - italic_ω ( italic_x ) ] italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) = 0 , (6b)

where we have defined,

μ=cM,𝜇𝑐𝑀\displaystyle\mu=cM,\hskip 5.0ptitalic_μ = italic_c italic_M , σ(x)=cS(r=cx),𝜎𝑥𝑐𝑆𝑟𝑐𝑥\displaystyle\sigma(x)=cS(r\!=\!cx),italic_σ ( italic_x ) = italic_c italic_S ( italic_r = italic_c italic_x ) , (7a)
ϵ=cE,italic-ϵ𝑐𝐸\displaystyle\epsilon=cE,\hskip 5.0ptitalic_ϵ = italic_c italic_E , ω(x)=cV(r=cx).𝜔𝑥𝑐𝑉𝑟𝑐𝑥\displaystyle\omega(x)=cV(r\!=\!cx).italic_ω ( italic_x ) = italic_c italic_V ( italic_r = italic_c italic_x ) . (7b)

Finally, we note that the Dirac orbitals satisfy the following normalization condition:

0(gnκ2(x)+fnκ2(x))𝑑x=1.superscriptsubscript0superscriptsubscript𝑔𝑛𝜅2𝑥superscriptsubscript𝑓𝑛𝜅2𝑥differential-d𝑥1\int_{0}^{\infty}\Big{(}g_{n\kappa}^{2}(x)+f_{n\kappa}^{2}(x)\Big{)}dx=1.∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_x ) + italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_x ) ) italic_d italic_x = 1 . (8)

II.2 Reduced Order Models

The main goal of this section is to illustrate how the use of a reduced order model helps generate efficient and accurate emulators for the solution of the Dirac equation. As mentioned in the Introduction, our aim is to speed up computations by approximating the solution of the high-fidelity problem by retaining only a handful of the most important components of the entire basis. But how does one identifies the most important components?

Our particular implementation of the reduced order model starts by obtaining exact solutions of the high-fidelity problem for a representative set of spherical nuclei. To facilitate the procedure, we approximate the self-consistent scalar and vector potentials generated from the covariant energy density functional FSUGarnet Chen and Piekarewicz (2015), by a 2-parameter Fermi (or Woods-Saxon) function, as illustrated in Fig.1 for 208Pb.

Refer to caption
Figure 1: Self-consistently generated scalar and time-like vector potentials for 208Pb, as predicted by the covariant energy density functional FSUGarnet Chen and Piekarewicz (2015). The dashed lines represent the corresponding 2-parameter Fermi fits that, after properly scaled, will be used to generate the reduced basis.
Nucleus Vssubscript𝑉sV_{\rm s}italic_V start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT Vvsubscript𝑉vV_{\rm v}italic_V start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT c𝑐citalic_c a𝑎aitalic_a λssubscript𝜆s\lambda_{\rm s}italic_λ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT λvsubscript𝜆v\lambda_{\rm v}italic_λ start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT β𝛽\betaitalic_β
O16superscriptO16{}^{16}\mathrm{O}start_FLOATSUPERSCRIPT 16 end_FLOATSUPERSCRIPT roman_O 442.8 360.6 2.605 0.475 5.862 4.747 5.495
Ca40superscriptCa40{}^{40}\mathrm{Ca}start_FLOATSUPERSCRIPT 40 end_FLOATSUPERSCRIPT roman_Ca 464.4 378.8 3.565 0.621 8.417 6.821 5.748
Ca48superscriptCa48{}^{48}\mathrm{Ca}start_FLOATSUPERSCRIPT 48 end_FLOATSUPERSCRIPT roman_Ca 438.4 358.4 3.948 0.511 8.774 7.169 7.739
Ni68superscriptNi68{}^{68}\mathrm{Ni}start_FLOATSUPERSCRIPT 68 end_FLOATSUPERSCRIPT roman_Ni 420.8 345.4 4.513 0.520 9.639 7.886 8.686
Zr90superscriptZr90{}^{90}\mathrm{Zr}start_FLOATSUPERSCRIPT 90 end_FLOATSUPERSCRIPT roman_Zr 411.5 335.7 5.034 0.507 10.52 8.550 9.935
Sn132superscriptSn132{}^{132}\mathrm{Sn}start_FLOATSUPERSCRIPT 132 end_FLOATSUPERSCRIPT roman_Sn 400.1 329.2 5.803 0.498 11.77 9.670 11.67
Pb208superscriptPb208{}^{208}\mathrm{Pb}start_FLOATSUPERSCRIPT 208 end_FLOATSUPERSCRIPT roman_Pb 402.4 330.2 6.795 0.565 13.86 11.357 12.02
Og302superscriptOg302{}^{302}\mathrm{Og}start_FLOATSUPERSCRIPT 302 end_FLOATSUPERSCRIPT roman_Og 381.0 313.5 7.845 0.553 15.17 12.447 14.19
Table 1: Optimal Woods-Saxon fit to the neutron scalar and vector potentials generated by the covariant energy density functional FSUGarnet Chen and Piekarewicz (2015). The strength of the scalar and vector potential is given in MeV, and the half-density radius and diffuseness in fm. The last three columns list the corresponding scaled, dimensionless parameters.

One then proceeds to scale both scalar and vector potentials as indicated in Eq.(7). That is,

S(r)=Vs1+e(rc)/aσ(x)=λs1+eβ(x1)𝑆𝑟subscript𝑉s1superscript𝑒𝑟𝑐𝑎𝜎𝑥subscript𝜆s1superscript𝑒𝛽𝑥1\displaystyle S(r)=\displaystyle{\frac{V_{\rm s}}{1+e^{(r-c)/a}}}% \longrightarrow\displaystyle{\sigma(x)=\frac{\lambda_{\rm s}}{1+e^{\,\beta(x-1% )}}}italic_S ( italic_r ) = divide start_ARG italic_V start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT ( italic_r - italic_c ) / italic_a end_POSTSUPERSCRIPT end_ARG ⟶ italic_σ ( italic_x ) = divide start_ARG italic_λ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT italic_β ( italic_x - 1 ) end_POSTSUPERSCRIPT end_ARG (9a)
V(r)=Vv1+e(rc)/aω(x)=λv1+eβ(x1),𝑉𝑟subscript𝑉v1superscript𝑒𝑟𝑐𝑎𝜔𝑥subscript𝜆v1superscript𝑒𝛽𝑥1\displaystyle V(r)=\displaystyle{\frac{V_{\rm v}}{1+e^{(r-c)/a}}}% \longrightarrow\displaystyle{\omega(x)=\frac{\lambda_{\rm v}}{1+e^{\,\beta(x-1% )}}},italic_V ( italic_r ) = divide start_ARG italic_V start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT ( italic_r - italic_c ) / italic_a end_POSTSUPERSCRIPT end_ARG ⟶ italic_ω ( italic_x ) = divide start_ARG italic_λ start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT italic_β ( italic_x - 1 ) end_POSTSUPERSCRIPT end_ARG , (9b)

where Vssubscript𝑉sV_{\rm s}italic_V start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT and Vvsubscript𝑉vV_{\rm v}italic_V start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT are the strengths of the scalar and vector potentials, respectively, c𝑐citalic_c is the half-density radius, and a𝑎aitalic_a is the diffuseness parameter. Once scaled, the resulting potentials depend linearly on the dimensionless strengths λssubscript𝜆s\lambda_{\rm s}italic_λ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT and λvsubscript𝜆v\lambda_{\rm v}italic_λ start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT, and non-linearly on the dimensionless ratio βc/a𝛽𝑐𝑎\beta\!\equiv\!c/aitalic_β ≡ italic_c / italic_a.

Refer to caption
Figure 2: (a) The upper component of the ground-state Dirac orbital generated by solving the dimensionful Dirac equation [Eq.(5)] for 48Ca, 90Zr, and 208Pb. (b) Same as (a) but now the Dirac orbitals are generated by solving the dimensionless Dirac equation [Eq.(6)]

Having defined the scaled potentials, one solves exactly the two-coupled differential equations given in Eqs.(6) by using a 4th-order Runge-Kutta algorithm for all nuclei listed in Table 1, with the exception of oganesson (Z=118𝑍118Z\!=\!118italic_Z = 118 and N=184𝑁184N\!=\!184italic_N = 184) that will be used later to test the efficacy of the emulator away from its training domain. As expected, the Dirac orbitals sharing the same quantum numbers are qualitatively similar. The simple rescaling implemented above makes the similarity between Dirac orbitals quantitative, a fact that is instrumental in generating a universal reduced basis. The impact of the scaling can be seen in Fig.2 where the left-hand panel shows the ground state 1s1/2(κ=1)1superscript𝑠12𝜅11s^{1/2}(\kappa\!=\!-1)1 italic_s start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT ( italic_κ = - 1 ) orbital generated by solving the original (dimensionful) Dirac equation. Instead, the right-hand panel shows the same orbitals after rescaling the Dirac equation to a dimensionless form. As much of the A𝐴Aitalic_A-dependence of the orbitals is mitigated by the scaling, the same reduced basis will be used to compute the single-energy spectrum of all spherical nuclei. Note that in the work reported in Ref. Giuliani et al. (2023), a reduced basis was created for each nucleus, and indeed, for each Dirac orbital. The availability of a universal reduced basis should further increase the efficiency of the emulator.

The last step in generating the (hopefully small) reduced basis invokes the Singular Value Decomposition (SVD) of an arbitrary matrix A𝐴Aitalic_A. As in our earlier work Anderson et al. (2022), we assemble all high fidelity solutions with the same value of κ𝜅\kappaitalic_κ into a rectangular A𝐴Aitalic_A matrix of dimension m×n𝑚𝑛m\times nitalic_m × italic_n and of rank r𝑟ritalic_r. That is,

A=UΣVT=u1σ1v1T+u2σ2v2T++urσrvrT𝐴𝑈Σsuperscript𝑉Tsubscript𝑢1subscript𝜎1superscriptsubscript𝑣1𝑇subscript𝑢2subscript𝜎2superscriptsubscript𝑣2𝑇subscript𝑢𝑟subscript𝜎𝑟superscriptsubscript𝑣𝑟𝑇A=U\Sigma V^{\rm T}=u_{1}\sigma_{1}v_{1}^{T}+u_{2}\sigma_{2}v_{2}^{T}+\ldots+u% _{r}\sigma_{r}v_{r}^{T}italic_A = italic_U roman_Σ italic_V start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT = italic_u start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT + italic_u start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT + … + italic_u start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT (10)

where U𝑈Uitalic_U and V𝑉Vitalic_V are m×m𝑚𝑚m\times mitalic_m × italic_m and n×n𝑛𝑛n\times nitalic_n × italic_n orthogonal matrices, respectively, and ΣΣ\Sigmaroman_Σ is an m×n𝑚𝑛m\times nitalic_m × italic_n “diagonal” matrix that encodes the rank of A𝐴Aitalic_A, namely, the dimensionality of the vector space. As such, the matrix A𝐴Aitalic_A can be decomposed as a sum of rank-1 matrices, as in Eq.(10), with the singular values σisubscript𝜎𝑖\sigma_{i}italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT arranged in descending order of importance. This criteria identifies the most important elements of the orthonormal reduced basis. Moreover, the singular values estimate the error likely to incur from truncations of the reduced basis. Finally, note that each vector that we input into the matrix is of dimension D=2Xmax𝐷2subscript𝑋maxD\!=\!2X_{\rm max}italic_D = 2 italic_X start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, where Xmaxsubscript𝑋maxX_{\rm max}italic_X start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT is the number of spatial grid points required to describe the Dirac orbitals: Xmaxsubscript𝑋maxX_{\rm max}italic_X start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT for the upper component and Xmaxsubscript𝑋maxX_{\rm max}italic_X start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT for the lower component.

II.3 The Reduced Basis Method

Having generated the reduced basis, we are now in a position to implement the reduced basis method approach to calculate the eigenstates of the Dirac Hamiltonian. The main virtue of such an approach is that whereas the solution of the high-fidelity model may be computationally demanding, the RBM transforms such challenging task into a diagonalization in a low-dimensional space spanned by the most important elements of the reduced basis. Having taken advantage of the spherical symmetry of the problem, the Dirac Hamiltonian matrix to be diagonalized—for each value of κ𝜅\kappaitalic_κ—is given by

H^(κ)=(μσ(x)+ω(x)ddx+κxddx+κxμ+σ(x)+ω(x)).^𝐻𝜅matrix𝜇𝜎𝑥𝜔𝑥𝑑𝑑𝑥𝜅𝑥𝑑𝑑𝑥𝜅𝑥𝜇𝜎𝑥𝜔𝑥\hat{H}(\kappa)=\begin{pmatrix}\mu-\sigma(x)+\omega(x)&\displaystyle{-\frac{d}% {dx}+\frac{\kappa}{x}}\\ \displaystyle{\frac{d}{dx}+\frac{\kappa}{x}}&-\mu+\sigma(x)+\omega(x)\end{% pmatrix}.over^ start_ARG italic_H end_ARG ( italic_κ ) = ( start_ARG start_ROW start_CELL italic_μ - italic_σ ( italic_x ) + italic_ω ( italic_x ) end_CELL start_CELL - divide start_ARG italic_d end_ARG start_ARG italic_d italic_x end_ARG + divide start_ARG italic_κ end_ARG start_ARG italic_x end_ARG end_CELL end_ROW start_ROW start_CELL divide start_ARG italic_d end_ARG start_ARG italic_d italic_x end_ARG + divide start_ARG italic_κ end_ARG start_ARG italic_x end_ARG end_CELL start_CELL - italic_μ + italic_σ ( italic_x ) + italic_ω ( italic_x ) end_CELL end_ROW end_ARG ) . (11)

In turn, matrix elements of the Dirac Hamiltonian are given by the following expression

m|H^(κ)|n=0(gmκ(x),fmκ(x))(μσ(x)+ω(x)ddx+κxddx+κxμ+σ(x)+ω(x))(gnκ(x)fnκ(x))𝑑x,quantum-operator-product𝑚^𝐻𝜅𝑛superscriptsubscript0subscript𝑔𝑚𝜅𝑥subscript𝑓𝑚𝜅𝑥matrix𝜇𝜎𝑥𝜔𝑥𝑑𝑑𝑥𝜅𝑥𝑑𝑑𝑥𝜅𝑥𝜇𝜎𝑥𝜔𝑥matrixsubscript𝑔𝑛𝜅𝑥subscript𝑓𝑛𝜅𝑥differential-d𝑥\langle m|\hat{H}(\kappa)|n\rangle={\int}_{0}^{\infty}\Big{(}g_{m\kappa}(x),f_% {m\kappa}(x)\Big{)}\begin{pmatrix}\mu-\sigma(x)+\omega(x)&\displaystyle{-\frac% {d}{dx}+\frac{\kappa}{x}}\\ \displaystyle{\frac{d}{dx}+\frac{\kappa}{x}}&-\mu+\sigma(x)+\omega(x)\end{% pmatrix}\begin{pmatrix}g_{n\kappa}(x)\vspace{3pt}\\ f_{n\kappa}(x)\end{pmatrix}dx,⟨ italic_m | over^ start_ARG italic_H end_ARG ( italic_κ ) | italic_n ⟩ = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( italic_g start_POSTSUBSCRIPT italic_m italic_κ end_POSTSUBSCRIPT ( italic_x ) , italic_f start_POSTSUBSCRIPT italic_m italic_κ end_POSTSUBSCRIPT ( italic_x ) ) ( start_ARG start_ROW start_CELL italic_μ - italic_σ ( italic_x ) + italic_ω ( italic_x ) end_CELL start_CELL - divide start_ARG italic_d end_ARG start_ARG italic_d italic_x end_ARG + divide start_ARG italic_κ end_ARG start_ARG italic_x end_ARG end_CELL end_ROW start_ROW start_CELL divide start_ARG italic_d end_ARG start_ARG italic_d italic_x end_ARG + divide start_ARG italic_κ end_ARG start_ARG italic_x end_ARG end_CELL start_CELL - italic_μ + italic_σ ( italic_x ) + italic_ω ( italic_x ) end_CELL end_ROW end_ARG ) ( start_ARG start_ROW start_CELL italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) end_CELL end_ROW start_ROW start_CELL italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) end_CELL end_ROW end_ARG ) italic_d italic_x , (12)

where (gnκ,fnκ)subscript𝑔𝑛𝜅subscript𝑓𝑛𝜅(g_{n\kappa},f_{n\kappa})( italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT , italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ) is the nthsubscript𝑛th{n\rm_{th}}italic_n start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT element of the reduced basis in the κ𝜅\kappaitalic_κ sector.

III Results

In this section we present in detail the path travelled to reach our ultimate goal, namely, the construction of a universal reduced basis for the efficient calibration of covariant energy density functionals. We decided to illustrate some of our missteps to underscore some of the subtleties behind the Dirac equation. Naturally, we started by assuming that the same approach employed in our previous publication to solve the Schrödinger equation Anderson et al. (2022) could also be used in the case of the Dirac equation.

III.1 The Naive Approach: Building a reduced basis
from only positive energy states

As in the case of the Schrödinger equation, the high fidelity model is solved exactly for a given set of model parameters and then one collects all bound orbitals with a given quantum number. In the present case, we use the bound neutron orbitals generated by all nuclei listed in Table 1, with the exception of Oganesson. This collection of non-orthogonal states is then filtered through an SVD routine to generate an optimal orthonormal basis arranged in order of importance. The reduced basis is then obtained by keeping only those orbitals with a condition number of at least 103superscript10310^{-3}10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT; that is, we keep orbitals whose singular value lies within 0.1% of the largest one. Note that the 103superscript10310^{-3}10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT value provides a rule-of-thumb estimate that worked well in our previous study Anderson et al. (2022).

Refer to caption
Figure 3: The bound single-particle spectrum of 90Zr. The black circles are the values obtained from the numerical Runge-Kutta integration, while the blue stars are the values obtained from the RBM. The dashed line indicates the approximate location of the Fermi energy; all bound states below the line are occupied where all the states above the line are empty.

Following such a procedure and using 90Zr as an example, we show in Fig.3 the entire single-neutron spectrum as generated by the high-fidelity model (black circles) as well as through the reduced basis method (blue stars). Note that the energies E𝐸Eitalic_E reported here include the neutron rest mass of M=940MeV𝑀940MeVM\!=\!940\,{\rm MeV}italic_M = 940 roman_MeV, so the corresponding binding energies are simply given by B=EM<0𝐵𝐸𝑀0B\!=\!E\!-\!M\!<\!0italic_B = italic_E - italic_M < 0. Although the diagonalization of the matrix within the reduced basis correctly computes all bound state energies—both occupied and unoccupied—many unphysical states also appear.

What may be the cause for the appearance of these unwanted “ghost” energies? The answer to this question lies in the nature of the Dirac spectrum. Indeed, the positive energy states of the Dirac equation—by themselves—do not form a complete basis. One of Dirac’s great insights was the discovery of negative-energy states that shortly after were identified as antiparticles. To shed light on how these ghost states emerge, we can regard the matrix diagonalization in the reduced basis as a variational approach that aims to select the N𝑁Nitalic_N optimal expansion parameters to describe the lowest energy state for each quantum number. That is, we can write

|ψ=n=1Nan|n,ket𝜓superscriptsubscript𝑛1𝑁subscript𝑎𝑛ket𝑛|\psi\rangle=\sum_{n=1}^{N}a_{n}|n\rangle,| italic_ψ ⟩ = ∑ start_POSTSUBSCRIPT italic_n = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT italic_a start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT | italic_n ⟩ , (13)

where |ψket𝜓|\psi\rangle| italic_ψ ⟩ is an eigenstate of the Dirac Hamiltonian, ansubscript𝑎𝑛a_{n}italic_a start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT are the expansion coefficients, and |nket𝑛|n\rangle| italic_n ⟩ is a member of the reduced basis introduced in Eq.(12). Although the structure of the positive energy states is accurately imprinted in the reduced basis (e.g., the upper components are significantly larger than the corresponding lower components) the Dirac Hamiltonian “knows” of the negative energy states. Hence, the extra ghost energies are an attempt to find the lowest energy solution. However, the attempt fails because bound states of negative energy can not be accurately expanded in terms of a reduced basis built from only positive energy states.

III.2 The Comprehensive Approach: Building a reduced basis
from positive and negative energy states

Faced with this dilemma, we now proceed to generate a new set of high-fidelity states, but this time by including both positive and negative energy states. As before, we use bound neutron orbitals generated by all nuclei listed in Table,̇1, with the exception of Oganesson. Although this involves a simple change in sign of the energy ϵitalic-ϵ\epsilonitalic_ϵ appearing in Eqs.(6), the impact of this change is dramatic. Whereas in the case of the positive energy states the large scalar and vector potentials combine to produce a fairly weak central potential, the result is entirely different for the negative energy states. In particular, the central component of the effective Schrödinger equation, given by

Uc(x)σ(x)+ϵμω(x),subscript𝑈c𝑥𝜎𝑥italic-ϵ𝜇𝜔𝑥U_{\rm c}(x)\approx-\sigma(x)+\frac{\epsilon}{\mu}\,\omega(x),italic_U start_POSTSUBSCRIPT roman_c end_POSTSUBSCRIPT ( italic_x ) ≈ - italic_σ ( italic_x ) + divide start_ARG italic_ϵ end_ARG start_ARG italic_μ end_ARG italic_ω ( italic_x ) , (14)

is fairly small for positive energy states because of the strong cancellation between large scalar and vector potentials. However, for negative energy states the strong cancellation is replaced by a strong enhancement, that leads to an explosion in the number of bound negative energy states. For example, in the case of 90Zr displayed in Fig.4, we find about 140 bound states—with only 15 of them having positive energy.

Refer to caption
Figure 4: The entire—both positive and negative—bound single-particle spectrum of 90Zr. The black circles are the values obtained from the numerical Runge-Kutta integration, while the blue stars are the values obtained from the RBM. The dashed line indicates the approximate location of the Fermi energy; all positive-energy bound states below the line are occupied whereas all the states above the line are empty.

While the inclusion of the entire spectrum of positive and negative energy states eliminated most of the unwanted ghost energies, some problems still remain. First, eliminating the ghost energies came at a very heavy price, as the number of basis states required to reach the necessary accuracy became enormous, reaching as high as 67 for some values of κ𝜅\kappaitalic_κ. Indeed, in order to eliminate all ghost energies, a very large number of basis states had to be retained, invalidating the idea of a “reduced” basis. Second, even at this heavy cost, we observe some residual ghost energies. Although most of the ghost energies that appeared in Fig.3 have been eliminated, a single d3/2superscript𝑑32d^{3/2}italic_d start_POSTSUPERSCRIPT 3 / 2 end_POSTSUPERSCRIPT ghost energy remains. Moreover, although most negative energy states are well reproduced, several unwanted ghost energies remain clearly visible. Note that no effort was made to explain the reason for the few remaining ghosts.

III.3 The Successful Approach: Building a reduced basis from
only positive energy states, but with a twist

In an effort to overcome the problem of relying on an overly large reduced basis which mitigates but does not entirely eliminate the appearance of ghost energies, we return to the original reduced basis built entirely from positive energy states. Although many unwanted ghost energies appear in Fig.3, the method seems to also reproduce the fifteen bound state energies obtained by solving the high-fidelity model. Hence, the task at hand consists on filtering out the ghost energies from the true energies. To do so, we implement the following two-fold prescription. First, we start by constructing a reduced basis from only positive energy state that we classify in terms of two quantum numbers (n,κ)𝑛𝜅(n,\kappa)( italic_n , italic_κ ), where n𝑛nitalic_n is the “principal” quantum number that accounts for the total number of nodes Giuliani et al. (2023). Second, once the Hamiltonian matrix is diagonalized within this reduced vector space, we only keep the eigenstate having the largest overlap with the most important element of the reduced basis. The most important element of the reduced basis, as quantified by its singular value, is guaranteed to have the same nodal structure as the exact eigenstate; the remaining elements of the reduced basis are simply in charge of performing the fine tuning.

Refer to caption
Figure 5: (a) Solutions of the high-fidelity model for the upper component of the 1d5/21superscript𝑑521d^{5/2}1 italic_d start_POSTSUPERSCRIPT 5 / 2 end_POSTSUPERSCRIPT orbital for 90Zr, 132Sn, and 208Pb. Also shown is the most important element of the reduced basis. (b) The resulting orthonormal basis obtained by filtering the high-fidelity solution through an SVD algorithm.

As an example on how we implement this procedure, we show in Fig. 5 how to generate the (n=1,κ=3)formulae-sequence𝑛1𝜅3(n\!=\!1,\kappa\!=\!-3)( italic_n = 1 , italic_κ = - 3 ) Dirac orbital, namely, the 1d5/21superscript𝑑521d^{5/2}1 italic_d start_POSTSUPERSCRIPT 5 / 2 end_POSTSUPERSCRIPT orbital containing one internal node. The left-hand panel in the figure displays high fidelity solutions for the upper (“large”) component of the 1d5/21superscript𝑑521d^{5/2}1 italic_d start_POSTSUPERSCRIPT 5 / 2 end_POSTSUPERSCRIPT orbital for 90Zr, 132Sn, and 208Pb. Note that this orbital is barely bound in 90Zr; see Fig.4. In turn, the right-hand panel in Fig.5 displays the optimal orthonormal reduced basis generated by SVD. As anticipated, the basis with the largest condition number has the same structure as the high-fidelity states. To underscore the similarity, we added this most important reduced-basis basis state to the left-hand panel.

Refer to caption
Figure 6: Upper (“large”) and lower (“small”) components of the 1d5/21superscript𝑑521d^{5/2}1 italic_d start_POSTSUPERSCRIPT 5 / 2 end_POSTSUPERSCRIPT Dirac orbital in 90Zr, as obtained from the numerical solution (RK) and from using a reduced basis (RB) method. The inset shows the projection of the Dirac orbital into the three reduced basis states.

Once the reduced basis has been obtained, one proceeds to diagonalize the “tiny” 3×3333\!\times\!33 × 3 Dirac Hamiltonian matrix. The result of the diagonalization procedure is displayed in Fig.6, which shows results for both the large upper g(x)𝑔𝑥g(x)italic_g ( italic_x ) and small lower f(x)𝑓𝑥f(x)italic_f ( italic_x ) components. The high fidelity solution obtained from using the Runge-Kutta algorithm is shown by the black solid line while the one obtained from the matrix diagonalization is depicted with the blue stars. Shown in the inset are the reduced-basis amplitudes as defined in Eq.(13). As expected, the first element of the reduced basis carries the largest weight, the second element provides most of the fine tuning, and the third element can be largely ignored. This procedure is applied to all nuclei in Table 1 and the results for four of them are displayed in Fig.7. Note that Og302superscriptOg302{}^{302}\,{\rm Og}start_FLOATSUPERSCRIPT 302 end_FLOATSUPERSCRIPT roman_Og was not included in the generation of the reduced basis, so its inclusion tests the performance of the model in extrapolating away from its training domain. Also note that some of the bound states in Og302superscriptOg302{}^{302}\,{\rm Og}start_FLOATSUPERSCRIPT 302 end_FLOATSUPERSCRIPT roman_Og do not contain an RB counterpart because of the obvious limitations of a reduced basis built only from lighter nuclei. Finally, the optimal reduced basis to be used in future calculations will benefit from the addition of Og302superscriptOg302{}^{302}\,{\rm Og}start_FLOATSUPERSCRIPT 302 end_FLOATSUPERSCRIPT roman_Og.

Refer to caption
Figure 7: The bound, positive-energy single-particle spectrum of 48Ca 90Zr, 208Pb, and 302Og. The black circles denote the values obtained from the numerical RK integration, while the blue stars are the values obtained from the RBM. The dashed line indicates the approximate location of the Fermi energy; all bound states below the line are occupied where all the states above the line are empty. The uncomputed energies in 302Og are a result of not having input functions to create a basis for those states.

We now summarize the advantages of implementing the two-fold procedure described above. By concentrating on orbitals classified according to the two quantum numbers (n,κ)𝑛𝜅(n,\kappa)( italic_n , italic_κ ), one can generate a small and efficient reduced basis for each channel. By virtue of the small dimensionality of the reduced basis, the diagonalization of the Dirac matrix is extremely fast. From the few eigenstates that emerge from the diagonalization, we are only interested in selecting one, namely, the one with the largest overlap with the most important element of the reduced basis. For all the bound orbitals below the Fermi surface that we have examined, the largest overlap exceeds 90%. Moreover, after rescaling the Dirac Hamiltonian according to Eq.(6), the reduced basis is “universal”, in the sense that it remains unchanged as one goes from O16superscriptO16{}^{16}\,{\rm O}start_FLOATSUPERSCRIPT 16 end_FLOATSUPERSCRIPT roman_O to Og302superscriptOg302{}^{302}\,{\rm Og}start_FLOATSUPERSCRIPT 302 end_FLOATSUPERSCRIPT roman_Og.

We conclude this section by making contact with our earlier publication Anderson et al. (2022). While the framework just described, namely, an emulator created with a distinct reduced basis for each set of (n,κ)𝑛𝜅(n,\kappa)( italic_n , italic_κ ) was successful, we realized that the filtering algorithm used in our earlier work consisting of a single basis for each κ𝜅\kappaitalic_κ can also work, as long as we ensure that the selected eigenvectors continue to have the largest overlap with the most relevant element of the reduced basis. Although the resulting Dirac matrix is larger, and thus takes longer to diagonalize, the advantage of creating only a single matrix for each κ𝜅\kappaitalic_κ results in an additional speedup factor of about 2.5, relative to the version using multiple smaller matrices for each κ𝜅\kappaitalic_κ.

III.4 Linearizing the Potential

The emulator that we have constructed so far speeds up the calculation of Dirac orbitals by nearly a factor of 150. While not insignificant, even such a considerable speedup factor often strains our computational resources. In particular, Bayesian calibration of covariant energy density functionals demand performing this sort of calculations thousands (if not millions) of times. As already mentioned, the diagonalization of the Dirac matrix in small vector spaces is very fast, so most of the computational time is spent computing the matrix elements given in Eq.(12). So if this bottleneck could be eliminated by precomputing all matrix elements, one would gain an additional advantage.

However, given that the scaled Woods-Saxon potentials are non-linear in the parameter β=c/a𝛽𝑐𝑎\beta\!=\!c/\!aitalic_β = italic_c / italic_a listed in Table 1, one would have to precompute matrix elements for each individual nucleus, thereby breaking the “universality” of the entire approach. Yet, one could mitigate this by linearizing the potential, thereby improving the dimensionality reduction of the problem. For example, one could write the dimensionless Woods-Saxon form given in Eq.(9) as follows:

u(x;β)=11+eβ(x1)q=1Qbq(β)uq(x),𝑢𝑥𝛽11superscript𝑒𝛽𝑥1superscriptsubscript𝑞1𝑄subscript𝑏𝑞𝛽subscript𝑢𝑞𝑥{u}(x;\beta)=\frac{1}{1+e^{\,\beta(x-1)}}\!\approx\!\sum_{q=1}^{Q}b_{q}(\beta)% {u}_{q}(x),italic_u ( italic_x ; italic_β ) = divide start_ARG 1 end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT italic_β ( italic_x - 1 ) end_POSTSUPERSCRIPT end_ARG ≈ ∑ start_POSTSUBSCRIPT italic_q = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_Q end_POSTSUPERSCRIPT italic_b start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_β ) italic_u start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_x ) , (15)

where the functions uq(x)subscript𝑢𝑞𝑥{u}_{q}(x)italic_u start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_x ) may be obtained following a procedure similar to the one adopted to construct the reduced basis. That is, one inputs into SVD a set of dimensionless Woods-Saxon forms spanning the range of values of β𝛽\betaitalic_β given in Table 1. Because SVD generates an optimal set of orthonormal functions, the nucleus-specific coefficients may be readily computed. That is,

bq(β)=0u(x;β)uq(x)𝑑x.subscript𝑏𝑞𝛽superscriptsubscript0𝑢𝑥𝛽subscript𝑢𝑞𝑥differential-d𝑥b_{q}(\beta)=\int_{0}^{\infty}{u}(x;\beta){u}_{q}(x)dx.italic_b start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_β ) = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_u ( italic_x ; italic_β ) italic_u start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_x ) italic_d italic_x . (16)

Hence, by linearizing the potential, the information specific to a given nucleus is encoded in the coefficients bq(β)subscript𝑏𝑞𝛽b_{q}(\beta)italic_b start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_β ), but matrix elements of uq(x)subscript𝑢𝑞𝑥{u}_{q}(x)italic_u start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_x ) between reduced basis states become independent of the nucleus under consideration. For example, matrix elements of the dimensionless vector potential ω(x;β)𝜔𝑥𝛽\omega(x;\beta)italic_ω ( italic_x ; italic_β ) may now be written according to Eq.(12) as follows:

mκ|ω^(β)|nκ=quantum-operator-product𝑚𝜅^𝜔𝛽𝑛𝜅absent\displaystyle\langle m\kappa|\hat{\omega}(\beta)|n\kappa\rangle=⟨ italic_m italic_κ | over^ start_ARG italic_ω end_ARG ( italic_β ) | italic_n italic_κ ⟩ = λv(β)mκ|u^(β)|nκsubscript𝜆v𝛽quantum-operator-product𝑚𝜅^𝑢𝛽𝑛𝜅\displaystyle\lambda_{\rm v}(\beta)\langle m\kappa|\hat{{u}}(\beta)|n\kappa\rangleitalic_λ start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT ( italic_β ) ⟨ italic_m italic_κ | over^ start_ARG italic_u end_ARG ( italic_β ) | italic_n italic_κ ⟩ (17)
=\displaystyle== λv(β)q=1Qbq(β)mκ|u^q|nκ.subscript𝜆v𝛽superscriptsubscript𝑞1𝑄subscript𝑏𝑞𝛽quantum-operator-product𝑚𝜅subscript^𝑢𝑞𝑛𝜅\displaystyle\lambda_{\rm v}(\beta)\sum_{q=1}^{Q}b_{q}(\beta)\langle m\kappa|% \hat{{u}}_{q}|n\kappa\rangle.italic_λ start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT ( italic_β ) ∑ start_POSTSUBSCRIPT italic_q = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_Q end_POSTSUPERSCRIPT italic_b start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_β ) ⟨ italic_m italic_κ | over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT | italic_n italic_κ ⟩ .

As such, the entire nuclear dependence is contained in the strength of the potential λv(β)subscript𝜆v𝛽\lambda_{\rm v}(\beta)italic_λ start_POSTSUBSCRIPT roman_v end_POSTSUBSCRIPT ( italic_β ) and in the bq(β)subscript𝑏𝑞𝛽b_{q}(\beta)italic_b start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_β ) coefficients. The remaining matrix elements, given by

mκ|u^q|nκ=0(gmκ(x)gnκ(x)+fmκ(x)fnκ(x))uq(x)𝑑x,quantum-operator-product𝑚𝜅subscript^𝑢𝑞𝑛𝜅superscriptsubscript0subscript𝑔𝑚𝜅𝑥subscript𝑔𝑛𝜅𝑥subscript𝑓𝑚𝜅𝑥subscript𝑓𝑛𝜅𝑥subscript𝑢𝑞𝑥differential-d𝑥\langle m\kappa|\hat{{u}}_{q}|n\kappa\rangle\!=\!\!\int_{0}^{\infty}\!\!\Big{(% }g_{m\kappa}(x)g_{n\kappa}(x)\!+\!f_{m\kappa}(x)f_{n\kappa}(x)\Big{)}{u}_{q}(x% )dx,⟨ italic_m italic_κ | over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT | italic_n italic_κ ⟩ = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( italic_g start_POSTSUBSCRIPT italic_m italic_κ end_POSTSUBSCRIPT ( italic_x ) italic_g start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) + italic_f start_POSTSUBSCRIPT italic_m italic_κ end_POSTSUBSCRIPT ( italic_x ) italic_f start_POSTSUBSCRIPT italic_n italic_κ end_POSTSUBSCRIPT ( italic_x ) ) italic_u start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT ( italic_x ) italic_d italic_x , (18)

are “universal”, so they can be precomputed offline and then stored for later use.

Once the potentials have been linearized and the relevant matrix elements computed and stored, all that remains is to build the matrix elements and diagonalize the various Dirac matrices in the reduced vector spaces. Using three expansion coefficients to linearize the potentials, the resulting emulator is about 12,000 times faster than the numerical RK integration method, and almost 100 times faster than using the reduced basis method in its original form. Although the linearization introduces a small additional error, the error is not visible when plotted against the exact numerical result. Hence the gains in computational expediency afforded by the linearization more than compensate for the slightly larger error.

IV Conclusions and Outlook

Emulators of high-fidelity models provide an efficient and accurate way to solve problems that often require long, computationally expensive numerical calculations. The Bayesian calibration of covariant energy density functionals falls into this category as it involves computing several nuclear properties for a variety of nuclei thousands of times. As part of the calibration protocol, one must solve the Dirac equation for all the occupied orbitals for all nuclei under consideration. In the high-fidelity model this is done by using the numerical RK integration method. Moreover, because the potentials generating the bound states depend themselves on the occupied orbitals, the problem becomes non-linear and must then be solved self-consistently, which exacerbates the computational demands.

In this paper we have extended our earlier work on the non-relativistic Schrödinger problem to the relativistic domain. As in the non-relativistic case, we built an efficient and accurate emulator by adopting a universal reduced basis capable of generating the entire single-particle spectrum of a variety of spherical nuclei without the need to redefine the basis. However, in contrast to the Schrödinger case, generating the relativistic reduced basis was complicated by the unavoidable existence of negative energy states. Although the problem was ultimately solved, we decided to document some of our tribulations, both because of the subtleties of the Dirac equation and to alert the reader of some our missteps. Ultimately, generating an efficient and robust basis required a few simple modifications relative to the strategy adopted in our earlier paper Anderson et al. (2022).

The successful reduced basis was built by generating high-fidelity solutions of the Dirac equations over a relative wide range of angular momentum κ𝜅\kappaitalic_κ. These solutions were then classified according to the quantum number κ𝜅\kappaitalic_κ. For each value of κ𝜅\kappaitalic_κ, an orthonormal reduced basis was generated using a singular value decomposition. Next, matrix elements of the Dirac Hamiltonian were generated for each κ𝜅\kappaitalic_κ and then the matrix was brought to diagonal form. From all the eigenvectors, one selects those having the largest overlap with one of the three most important reduced basis states, namely, with the ones having the largest singular values. This guarantees that all unwanted “ghosts” states that contaminated the spectrum disappear. We note that for all Dirac orbitals and for all nuclei, accurate results were obtained by diagonalizing at most a 12×12121212\!\times\!1212 × 12 matrix. Moreover, in most of the cases, the most important element of the reduced basis carried at least 90% of the weight. With this implementation, we saw speed ups of about a factor of 150, relative to the high-fidelity model.

Whereas a factor of 150 is significant, we identified the evaluation of the matrix elements as the most serious bottleneck in the problem. This bottleneck develops because the scalar and vector potentials depend non-linearly on the parameter β𝛽\betaitalic_β, which is different for each nucleus. We mitigated the problem by invoking a linearization scheme. The linearization improves greatly the dimensionality reduction of the problem by expanding the non-linear potential in a set of basis states that are independent of the nucleus and a corresponding set of coefficients that carry the nuclear information. In this manner, matrix elements of the nucleus-independent basis states may be precomputed and stored, leading to an even greater computational expediency. Indeed, precomputing matrix elements of the potential in the small reduced vector spaces yields an additional speed-up factor of 100 relative to the RBM in its original form. That is, by implementing these methods the calculations can be performed thousands of times faster, without sacrificing accuracy.

In the future, we aim to mitigate the last remaining bottleneck associated with the iterative procedure required to solve the mean-field problem self-consistently. Indeed, besides a well-motivated reduced basis that captures the essential physics, one also requires an efficient framework to solve the underlying set of dynamical equations. For a system of linear differential equation as the one considered here, the solution simply involves direct diagonalization of the Hamiltonian matrix in the reduced-basis space. However, for the self-consistent solution of the mean-field problem, one needs to solve a non-linear set of differential equations. In this case, rather than solving the self-consistent problem by iteration, one invokes the Galerkin projection approach, an enormously efficient framework that involves solving once a non-linear set of algebraic equations Quarteroni et al. (2015). Although without the benefit of a universal reduced basis, significant progress along this direction has already been accomplished in Ref. Giuliani et al. (2023).

Acknowledgements.
We are grateful to Pablo Giuliani for many useful and insightful discussions. This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Nuclear Physics under Award Number DE-FG02-92ER40750 and under the STREAMLINE collaboration award DE-SC0024646.

References

  • PRA-Editors (2011) PRA-Editors, Phys. Rev. A 83, 040001 (2011).
  • Reinhard and Nazarewicz (2010) P.-G. Reinhard and W. Nazarewicz, Phys. Rev. C81, 051303(R) (2010).
  • Kortelainen et al. (2010) M. Kortelainen, T. Lesinski, J. More, W. Nazarewicz, J. Sarich, et al., Phys. Rev. C82, 024313 (2010).
  • Fattoyev and Piekarewicz (2011) F. Fattoyev and J. Piekarewicz, Phys. Rev. C84, 064302 (2011).
  • Reinhard et al. (2013) P.-G. Reinhard, J. Piekarewicz, W. Nazarewicz, B. Agrawal, N. Paar, et al., Phys.Rev. C88, 034325 (2013).
  • Dobaczewski et al. (2014) J. Dobaczewski, W. Nazarewicz, and P.-G. Reinhard, J. Phys. G41, 074001 (2014).
  • Piekarewicz et al. (2015) J. Piekarewicz, W.-C. Chen, and F. Fattoyev, J. Phys. G42, 034018 (2015).
  • Chen and Piekarewicz (2014) W.-C. Chen and J. Piekarewicz, Phys. Rev. C90, 044305 (2014).
  • Wesolowski et al. (2016) S. Wesolowski, N. Klco, R. Furnstahl, D. Phillips, and A. Thapaliya, J. Phys. G 43, 074001 (2016).
  • Yuksel et al. (2019) E. Yuksel, T. Marketin, and N. Paar, Phys. Rev. C99, 034318 (2019).
  • Drischler et al. (2020a) C. Drischler, R. Furnstahl, J. Melendez, and D. Phillips, Phys. Rev. Lett. 125, 202702 (2020a).
  • Drischler et al. (2020b) C. Drischler, J. A. Melendez, R. J. Furnstahl, and D. R. Phillips, Phys. Rev. C 102, 054315 (2020b).
  • Furnstahl et al. (2020) R. J. Furnstahl, A. J. Garcia, P. J. Millican, and X. Zhang, Phys. Lett. B 809, 135719 (2020).
  • Athanassopoulos et al. (2004) S. Athanassopoulos, E. Mavrommatis, K. A. Gernoth, and J. W. Clark, Nucl. Phys. A743, 222 (2004).
  • Akkoyun et al. (2013) S. Akkoyun, T. Bayram, S. O. Kara, and A. Sinan, J. Phys. G40, 055106 (2013).
  • Bayram et al. (2014) T. Bayram, S. Akkoyun, and S. O. Kara, Annals of Nuclear Energy 63, 172 (2014).
  • Gernoth et al. (1993) K. Gernoth, J. Clark, J. Prater, and H. Bohr, Phys. Lett. B 300, 1 (1993).
  • Utama et al. (2016a) R. Utama, J. Piekarewicz, and H. B. Prosper, Phys. Rev. C93, 014311 (2016a).
  • Utama et al. (2016b) R. Utama, W.-C. Chen, and J. Piekarewicz, J. Phys. G, 014311 (2016b).
  • Utama and Piekarewicz (2018) R. Utama and J. Piekarewicz, Phys. Rev. C97, 014306 (2018).
  • Neufcourt et al. (2019) L. Neufcourt, Y. Cao, W. Nazarewicz, E. Olsen, and F. Viens, Phys. Rev. Lett. 122, 062502 (2019).
  • Neufcourt et al. (2020) L. Neufcourt, Y. Cao, S. A. Giuliani, W. Nazarewicz, E. Olsen, and O. B. Tarasov, Phys. Rev. C 101, 044307 (2020).
  • Lovell et al. (2022) A. E. Lovell, A. T. Mohan, T. M. Sprouse, and M. R. Mumpower, Phys. Rev. C 106, 014305 (2022).
  • Saito et al. (2024) Y. Saito, I. Dillmann, R. Kruecken, M. R. Mumpower, and R. Surman, Phys. Rev. C 109, 054301 (2024).
  • Frame et al. (2018) D. Frame, R. He, I. Ipsen, D. Lee, D. Lee, and E. Rrapaj, Phys. Rev. Lett. 121, 032501 (2018).
  • König et al. (2020) S. König, A. Ekström, K. Hebeler, D. Lee, and A. Schwenk, Phys. Lett. B 810, 135814 (2020).
  • Drischler et al. (2021) C. Drischler, M. Quinonez, P. G. Giuliani, A. E. Lovell, and F. M. Nunes, Phys. Lett. B 823, 136777 (2021).
  • Bonilla et al. (2022) E. Bonilla, P. Giuliani, K. Godbey, and D. Lee, Phys. Rev. C 106, 054322 (2022).
  • Melendez et al. (2022) J. A. Melendez, C. Drischler, R. J. Furnstahl, A. J. Garcia, and X. Zhang, J. Phys. G 49, 102001 (2022).
  • Giuliani et al. (2023) P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Front. Phys. 10, 1054524 (2023).
  • Anderson et al. (2022) A. L. Anderson, G. L. O’Donnell, and J. Piekarewicz, Phys. Rev. C 106, L031302 (2022).
  • Odell et al. (2024) D. Odell, P. Giuliani, K. Beyer, M. Catacora-Rios, M. Y.-H. Chan, E. Bonilla, R. J. Furnstahl, K. Godbey, and F. M. Nunes, Phys. Rev. C 109, 044612 (2024).
  • (33) K. Godbey, P. Giuliani, E. Bonilla, E. Flynn, D. Odell, K. Beyer, D. Lay, D. Figueroa, R. Garg, and M. Campbell, Dimensionality reduction in nuclear physics, https://dr.ascsn.net/.
  • Brunton and Kutz (2019) S. L. Brunton and J. N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control (Cambridge University Press, 2019).
  • Benner et al. (2020) P. Benner, W. Schilders, S. Grivet-Talocia, A. Quarteroni, G. Rozza, and L. Silveira, Model Order Reduction: Applications (De Gruyter, 2020), vol. 3.
  • Quarteroni et al. (2015) A. Quarteroni, A. Manzoni, and F. Negri, Reduced Basis Methods for Partial Differential Equations: An Introduction, vol. 92 (Springer, 2015).
  • Hesthaven et al. (2016) J. S. Hesthaven, G. Rozza, and B. Stamm, Certified Reduced Basis Methods for Parametrized Partial Differential Equations (Springer, Switzerland, 2016), vol. 590.
  • Bohr and Mottelson (1998) A. Bohr and B. R. Mottelson, Nuclear Structure (World Scientific Publishing Company, New Jersey, 1998), vol. I: Single-particle motion.
  • Todd-Rutel and Piekarewicz (2005) B. G. Todd-Rutel and J. Piekarewicz, Phys. Rev. Lett 95, 122501 (2005).
  • Fattoyev et al. (2010) F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and G. Shen, Phys. Rev. C82, 055803 (2010).
  • Chen and Piekarewicz (2015) W.-C. Chen and J. Piekarewicz, Phys. Lett. B748, 284 (2015).
  • Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  • Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  • Dirac (1982) P. Dirac, The Principles of Quantum Mechanics (Clarendon Press, Oxford, 1982).
  • Horowitz and Serot (1981) C. J. Horowitz and B. D. Serot, Nucl. Phys. A368, 503 (1981).
  • Serot and Walecka (1986) B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. 16, 1 (1986).
  • Sakurai (1967) J. J. Sakurai, Advanced Quantum Mechanics (Pearson Education, 1967).