A universal reduced basis for the calibration of covariant energy density functionals
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:
(1) |
where and are Dirac matrices, is the momentum operator, and is the bare nucleon mass. In turn, and 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 and its projection Sakurai (1967). That is,
(2) |
where is the principal quantum number and the spin spherical harmonics are obtained by coupling the orbital angular momentum to the intrinsic nucleon spin-1/2 to obtain a total angular momentum . That is,
(3) |
In particular, the orbital angular momentum and the total angular momentum are obtained from the generalized angular momentum as follows:
(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:
(5a) | |||
(5b) |
In search of a universal basis, it is convenient to rescale the above equations by an intrinsic distance scale characteristic of each individual nucleus. Hence, by defining the dimensionless distance as , the coupled Dirac equation may now be written as
(6a) | |||
(6b) |
where we have defined,
(7a) | |||||
(7b) |
Finally, we note that the Dirac orbitals satisfy the following normalization condition:
(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.
Nucleus | |||||||
---|---|---|---|---|---|---|---|
442.8 | 360.6 | 2.605 | 0.475 | 5.862 | 4.747 | 5.495 | |
464.4 | 378.8 | 3.565 | 0.621 | 8.417 | 6.821 | 5.748 | |
438.4 | 358.4 | 3.948 | 0.511 | 8.774 | 7.169 | 7.739 | |
420.8 | 345.4 | 4.513 | 0.520 | 9.639 | 7.886 | 8.686 | |
411.5 | 335.7 | 5.034 | 0.507 | 10.52 | 8.550 | 9.935 | |
400.1 | 329.2 | 5.803 | 0.498 | 11.77 | 9.670 | 11.67 | |
402.4 | 330.2 | 6.795 | 0.565 | 13.86 | 11.357 | 12.02 | |
381.0 | 313.5 | 7.845 | 0.553 | 15.17 | 12.447 | 14.19 |
One then proceeds to scale both scalar and vector potentials as indicated in Eq.(7). That is,
(9a) | |||
(9b) |
where and are the strengths of the scalar and vector potentials, respectively, is the half-density radius, and is the diffuseness parameter. Once scaled, the resulting potentials depend linearly on the dimensionless strengths and , and non-linearly on the dimensionless ratio .
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 ( and ) 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 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 -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 . As in our earlier work Anderson et al. (2022), we assemble all high fidelity solutions with the same value of into a rectangular matrix of dimension and of rank . That is,
(10) |
where and are and orthogonal matrices, respectively, and is an “diagonal” matrix that encodes the rank of , namely, the dimensionality of the vector space. As such, the matrix can be decomposed as a sum of rank-1 matrices, as in Eq.(10), with the singular values 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 , where is the number of spatial grid points required to describe the Dirac orbitals: for the upper component and 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 —is given by
(11) |
In turn, matrix elements of the Dirac Hamiltonian are given by the following expression
(12) |
where is the element of the reduced basis in the 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 ; that is, we keep orbitals whose singular value lies within 0.1% of the largest one. Note that the value provides a rule-of-thumb estimate that worked well in our previous study Anderson et al. (2022).
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 reported here include the neutron rest mass of , so the corresponding binding energies are simply given by . 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 optimal expansion parameters to describe the lowest energy state for each quantum number. That is, we can write
(13) |
where is an eigenstate of the Dirac Hamiltonian, are the expansion coefficients, and 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 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
(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.
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 . 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 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 , where 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.
As an example on how we implement this procedure, we show in Fig. 5 how to generate the Dirac orbital, namely, the orbital containing one internal node. The left-hand panel in the figure displays high fidelity solutions for the upper (“large”) component of the 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.
Once the reduced basis has been obtained, one proceeds to diagonalize the “tiny” Dirac Hamiltonian matrix. The result of the diagonalization procedure is displayed in Fig.6, which shows results for both the large upper and small lower 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 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 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 .
We now summarize the advantages of implementing the two-fold procedure described above. By concentrating on orbitals classified according to the two quantum numbers , 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 to .
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 was successful, we realized that the filtering algorithm used in our earlier work consisting of a single basis for each 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 results in an additional speedup factor of about 2.5, relative to the version using multiple smaller matrices for each .
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 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:
(15) |
where the functions 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 given in Table 1. Because SVD generates an optimal set of orthonormal functions, the nucleus-specific coefficients may be readily computed. That is,
(16) |
Hence, by linearizing the potential, the information specific to a given nucleus is encoded in the coefficients , but matrix elements of between reduced basis states become independent of the nucleus under consideration. For example, matrix elements of the dimensionless vector potential may now be written according to Eq.(12) as follows:
(17) | ||||
As such, the entire nuclear dependence is contained in the strength of the potential and in the coefficients. The remaining matrix elements, given by
(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 . These solutions were then classified according to the quantum number . For each value of , an orthonormal reduced basis was generated using a singular value decomposition. Next, matrix elements of the Dirac Hamiltonian were generated for each 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 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 , 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).