Microscopically based energy density functionals for nuclei using the density matrix expansion

TL;DR

This paper develops a microscopically motivated energy density functional for nuclei using the density matrix expansion, aiming to improve accuracy and computational efficiency in nuclear structure calculations.

Contribution

It introduces an approximation of Yukawa functions with Gaussian sums to reduce computational cost in density functional calculations based on chiral nuclear forces.

Findings

Gaussian approximation achieves acceptable accuracy with fewer functions.

The method reduces computational memory requirements.

The approach enhances the practicality of microscopic energy density functionals.

Abstract

While ab initio many-body techniques have been able to successfully describe the properties of light and intermediate mass nuclei based on chiral effective field theory interactions, neutron-rich nuclei still remain out of reach for these methods. Conversely, energy density functional approaches can be used to calculate properties of heavy nuclei but rely mostly on phenomenological interactions. A usable form of the nuclear energy density functional that is rooted in the modern theory of nuclear forces was presented recently. The first component of this new set of functionals corresponds to the direct part (Hartree term) of the expectation value of local chiral potentials on a Slater determinant. The exchange term, which is a functional of the non-local density, is transformed into a local functional by applying the density matrix expansion. In order to reduce the computational cost due to the direct implementation of non-separable, local interactions in the Hartree term, we use an approximation to represent the regularized Yukawa functions in terms of a sum of (separable) Gaussian functions. These proceedings analyze the accuracy of such an approximation in terms of the number of Gaussian functions and look for an optimal value that gives an acceptable level of accuracy while maintaining the computational memory requirements in a many-body calculation as low as possible.

Figures

Tables

Table 1. Table 1: Root mean square deviation between the two central terms V C subscript 𝑉 𝐶 V_{C} and W C subscript 𝑊 𝐶 W_{C} described in  [ 17 ] and the Gaussian expansion in equation ( 6 ) for different number of Gaussian functions.

\br$N$	$V_c$	$W_c$
	MeV	MeV
\mr3	0.0229	0.0828
4	0.0083	0.0155
5	0.0013	0.0039
6	0.0013	0.0018
7	0.0010	0.0017
8	0.0010	0.0012
\br

Equations

$$E^{\mathrm{NN}}[\rho ]=\frac{1}{2}\sum _{ij}⟨ij\left|V^{\mathrm{NN}}(1-P_{\sigma }{P}_{\tau }{P}_r)\right|kl⟩{\rho }_{ki}{\rho }_{lj},$$

$$E_{\mathrm{Hartree}}^{\mathrm{NN}}$$

$$E_{\mathrm{Fock}}^{\mathrm{NN}}$$

$$\rho \left(\mathbf{R}+\frac{\mathbf{r}}{2},\mathbf{R}-\frac{\mathbf{r}}{2}\right)\approx {\Pi }_0^{\rho }(k_{F}r)\rho (\mathbf{R})+\frac{r^2}{6}{\Pi }_2^{\rho }(k_{F}r)\left[\frac{1}{4}\Delta \rho (\mathbf{R})-\tau (\mathbf{R})+\frac{3}{5}{k}_F^2\rho (\mathbf{R})\right],$$

$$\tau (\mathbf{r})=\nabla \cdot {\nabla }^{\prime }\rho (\mathbf{r},{\mathbf{r}}^{\prime }){|}_{\mathbf{r}={\mathbf{r}}^{\prime }}.$$

$$V_{\chi }(r)\to \tilde{V}(r)=\sum _{i=1}^N{V}_i{e}^{-r^2/{\mu }_i}$$

$$f(r)={\left[1-\exp \left(-\frac{r^2}{R_c^2}\right)\right]}^n.$$

$$V_N=-\sum _{i=1}^{N-1}{V}_i.$$

Full text

\setabbreviationstyle

[acronym]long-short \glssetcategoryattributeacronymnohyperfirsttrue

LLNL-PROC-769943

1 Introduction

One of the overarching goals of nuclear physics is to find a theory capable of describing nuclear structure and reactions in terms of the underlying interactions between quarks and gluons. The description of the interaction between nucleons in \glsxtrprotectlinksChiral Effective Field Theory ($\chi$-EFT) provides the desired microscopic underpinning and also exhibits systematic, order-by-order, improvement [1, 2, 3, 4]. Ab initio methods employ these realistic interactions and calculate properties of nuclear systems by directly considering the interaction between each pair of nucleons. This approach has been successful in describing light and medium nuclei (up to $Z=50$) [5]. The same systematic, order-by-order, improvement from the interaction can be observed in the description of nuclear properties calculated with ab initio methods [6]. This allows keeping track of the predictive power of these type of theoretical calculations when looking at properties that have not been measured yet. However heavy systems, particularly neutron-rich systems, remain out of reach for ab initio methods. These neutron-rich systems are of great interest in view of upcoming measurements in flagship facilities like \glsxtrprotectlinksthe Facility for Rare Isotope Beams (FRIB) and for the description of astrophysical phenomena like the r-process [7, 8].

In contrast to ab initio methods, \glsxtrprotectlinksEnergy Density Functional (EDF) methods, also known as self-consistent mean-field approaches, substitute the interaction due to all pairs of nucleons by an “average” of the interaction with all the other nucleons [9, 10]. This substitution significantly reduces the computational cost of many-body nuclear calculations and allows estimating properties of heavy systems such as the neutron-rich nuclei mentioned earlier. However, the form of the interaction used in \glsxtrprotectlinksEDF approaches is usually constructed ad hoc in order to simplify calculations, and parameters are adjusted to reproduce experimental nuclear properties such as masses, charge radii or spin-orbit splittings. This leads to effective interactions of rather phenomenological nature whose predictive power, when extrapolated to properties beyond the ones used to determine its parameters, can be hard to maintain or even estimate.

Following earlier studies [11], we have recently proposed a new method to determine an \glsxtrprotectlinksEDF from microscopic nuclear potentials derived from \glsxtrprotectlinks$\chi$-EFT through a tool known as the \glsxtrprotectlinksDensity Matrix Expansion (DME). In practice, the new \glsxtrprotectlinksEDF looks like a phenomenological Skyrme-like contact interaction to which is added a more microscopic component which takes the form of density-dependent couplings resulting from applying the \glsxtrprotectlinksDME to the non-local one-body density [12]. This approach encodes the long-range physics contained in the \glsxtrprotectlinks$\chi$-EFT potential. The phenomenological Skyrme-like part accounts for the short-range physics and its parameters need to be readjusted in order to reproduce experimental nuclear structure data.

The computational implementation of these semi-phenomenological interactions involves using certain numerical approximations (beyond the \glsxtrprotectlinksDME itself). In particular, the contribution of the long-range part of the chiral potential to the nuclear binding energy can be efficiently computed by approximating the Yukawa potentials with a sum of Gaussian functions. This approximation allows to recast the \glsxtrprotectlinks$\chi$-EFT-derived potential in a form analogous to the purely phenomenological Gogny interaction, taking advantage of fact that this type of interaction is already implemented in several many-body \glsxtrprotectlinksDensity Functional Theory (DFT) codes. Of course the implementation of any approximation requires to ensure a desired level of accuracy in order to avoid the introduction and propagation of numerical errors in all kinds of calculations. These proceedings analyze the convergence of the sum of Gaussian functions approximation to the central terms in the chiral potential at \glsxtrprotectlinksnext to leading order (NLO) and \glsxtrprotectlinksnext to next to leading order (N2LO). Section 2 gives a brief overview of the implementation of semi-phenomenological \glsxtrprotectlinksEDF introduced in [12]. Following that, section 3 presents an analysis of the accuracy of the Gaussian expansion depending on the number of terms included. Conclusions and outlook are discussed in section 4.

2 Semi-phenomenological EDF

Our many-body theoretical framework is the \glsxtrprotectlinksHartree-Fock-Bogoliubov (HFB) formalism. In this framework, the many-body wave function is taken as a quasiparticle vacuum [9]. The total energy is expressed as a functional of the generalized density, or equivalently, as a functional of the one-body density matrix $\rho$ and the pairing tensor $\kappa$. The full \glsxtrprotectlinksEDF can often be broken down into a sum of several terms, one that depends exclusively on $\rho$ (the particle-hole channel) and another one that also involves $\kappa$ and is only active when pairing correlations are present (the particle-particle channel).

The \glsxtrprotectlinksEDF in the particle-hole channel is often constructed by representing the many-body wave function by a Slater determinant and using the Wick theorem to express the energy due to many-body interactions as a functional of the density of particles in the system. In the simplest case of two-body interactions, this leads to writing

[TABLE]

with $\rho_{ij}$ the matrix elements of the one-body density matrix on an arbitrary basis of the single-particle Hilbert space, $P_{\sigma}$, $P_{\tau}$ and $P_{r}$ the usual spin-, isospin- and position-exchange operators and $V^{\rm NN}$ the two-body nuclear potential. The anti-symmetrization operator $(1-P_{\sigma}P_{\tau}P_{r})$ results in direct and exchange contributions, also known as the Hartree and Fock energies. After transforming the spatial coordinates to relative $(\bm{r})$ and center-of-mass $(\bm{R})$ coordinates, the Hartree and Fock energies become

[TABLE]

where the potential bracket includes spin and isospin states, the traces refer to summations over those quantum numbers and the local density matrix is $\rho(\bm{x})=\rho(\bm{x},\bm{x})$; for details on how to incorporate three-body forces into this formalism, see refs. [13, 12]. In a purely phenomenological approach the spatial form of the two-body interaction $V^{\rm NN}$ is partly chosen to simplify the calculation of these energy terms. Popular choices include delta functions for Skyrme functionals [14, 15] and Gaussian functions for the Gogny interaction [16]. Note that by representing the many-body wave function with a Slater determinant, the \glsxtrprotectlinksEDF approach assumes independent particles. Many-body correlations are recovered by a combination of two mechanisms: the spontaneous breaking of symmetries of the nuclear potential and the calibration of the parameters of the energy-density dependent functional to experimental data [9].

In order to anchor the \glsxtrprotectlinksEDF approach to a more microscopic theory and a systematic, order-by-order, description of nuclear forces, the new microscopic functionals presented in [12] start from the local chiral potentials up to \glsxtrprotectlinksN2LO with and without $\Delta$ excitations as outlined in [17, 18]. A direct substitution of these potentials in equations (2) and (3) instead of the purely phenomenological Gogny or Skyrme interactions results in computationally more expensive numerical integrals that would need to be calculated at every \glsxtrprotectlinksHFB iteration. To sidestep the increase in computational cost two approximations are made. The spatial form factor of the potential in equation (2) is approximated by a sum of Gaussian functions whose strength and range parameters are adjusted to reproduce the original chiral interaction. The use of Gaussian functions allows treating the Hartree term in the same fashion as a Gogny interaction and make use of previously benchmarked many-body \glsxtrprotectlinksDFT codes in which Gogny interactions have already been implemented, such as HFBTHOv3 [19]. The second approximation is to employ the \glsxtrprotectlinksDME formalism to calculate the Fock energy term. The \glsxtrprotectlinksDME expands the non-local density in equation (3) in terms of a local density and derivatives of that local density as

[TABLE]

where the $\Pi$ functions are proportional to a reduced spherical Bessel function, $k_{F}$ is the usual Fermi momentum and $\tau(\bm{r})$ is the kinetic density given by

[TABLE]

As a result, the \glsxtrprotectlinksDME turns the finite-range part of the local, \glsxtrprotectlinks$\chi$-EFT chiral potential into a Skyrme-like functional, the main difference with traditional Skyrme functionals being that all coupling constants are now functionals of the local density.

The procedure briefly outlined above results in a semi-phenomenological approach in which the re-parametrization of the Skyrme-like part encodes the zero-range, i.e., high-energy, part of the chiral potential and represents the short-range interaction, while the finite-range part of the interaction provides the order-by-order improvement to the description of many-body potentials. With this approach it is in principle possible to construct a scalable framework capable of describing neutron-rich nuclei with a series of systematic improvements.

3 Numerical accuracy of the approximation for the Hartree term

As mentioned in the previous section, the implementation of the microscopically constrained density functionals requires approximating the chiral potential in equation (2) as a sum of Gaussian functions,

[TABLE]

It is crucial then to ensure that this Gaussian expansion gives a correct description of the actual chiral interaction in terms of Yukawa potentials.

An important feature of the local chiral potentials of [17, 18] is that the potential contains a regulator in order to avoid divergences at short distances. The regulator has the form

[TABLE]

The potentials implemented in [12] in particular use $R_{c}=1.0$ fm and $n=6$. This regulator ensures that the interaction vanishes as $r\rightarrow 0$. A simple way to reproduce this kind of behavior in the sum of Gaussian functions in equation (6) is to impose the condition,

[TABLE]

The choice for the number of Gaussian functions used in the approximation is dictated by the necessary compromise between precision and available resources. Using too few Gaussian functions causes a poor approximation that will introduce numerical errors in the \glsxtrprotectlinksHFB many-body calculation; using too many Gaussian functions significantly increases the necessary memory to store all non-zero one-body density matrix elements with the possibility of going beyond the amount of memory available in current high performance computing architectures. With these restrictions in mind it is worth looking at the difference between the actual chiral interaction and equation (6) as a function of the number of Gaussian terms included in the approximation. While these proceedings focus in particular on the $W_{C}^{\rm NLO}$ and $V_{C}^{\rm N2LO}$ terms in the chiral interaction outlined in [17], the analysis presented here applies as well to the other central terms at different orders with and without the inclusion of $\Delta$ excitations.

For each of the two interactions mentioned above, we study six different approximations that go from using only three Gaussian functions up to a total of eight. In each approximation the $V_{i}$ and $\mu_{i}$ parameters are adjusted, imposing the constraint of equation (8), to reproduce the chiral interaction as well as possible. While a regularized approximation ($\tilde{V}\rightarrow 0$ as $r\rightarrow 0$) can already be constructed using only two Gaussian functions, a reasonable approximation is obtained starting at $N=3$. The same type of analysis for the central isospin term $W_{C}^{\rm NLO}$ is presented in figure 2. Although the curves in the right panel show a reduction in the difference between the chiral potential and the sum of Gaussian functions approximation from $N=5$ to $N=6$, the order for magnitude remains the same all the way up to $N=8$. Considering that the computational memory requirements would be considerably increased by an extra Gaussian function it is preferred to keep the approximation at $N=5$ as the optimal compromise between numerical accuracy and memory requirements.

Although the figures 1 and 2 give visual evidence that increasing the number of Gaussian functions beyond $N=5$ adds little to no improvement, it is preferred to quantify this statement via a numerical calculation. To estimate the average discrepancy between the chiral interaction and the Gaussian approximation (6), we calculate the root mean square deviation using a fine grid sampling in radius every $0.01$ fm from $r=0$ fm to $r=10$ fm. The numerical results are shown in table 1 and quantifies what was already suggested in the previous figures; there is little to no change in the error beyond $N=5$.

4 Conclusions and outlook

The implementation of microscopically-constrained \glsxtrprotectlinksEDF allows performing many-body calculations of large neutron-rich systems with a systematic order-by-order improvement determined by \glsxtrprotectlinks$\chi$-EFT. The performance of this new family of functionals has shown an improved description of nuclear masses when compared to other state-of-the-art functionals like UNEDF2 [12]. This milestone paves the way towards performing calculations relevant to astrophysics, such as simulations of the r-process, with nuclear models with a firm microscopic underpinning stemming from quarks and gluons. Given the potential of these new functionals, it is important to verify that the approximations employed for the implementation have a reasonable level of accuracy in order to avoid compromising the numerical precision of upcoming calculations. These proceedings focused on the numerical accuracy of the sum of Gaussian functions used to represent the central terms of the chiral potential in the calculation of the Hartree energy term. While an infinite number of Gaussian functions could potentially guarantee numerically exact representation of the regularized Yukawa potentials in the chiral interaction, computational resources are limited and a large number of Gaussian functions results in a prohibitively large amount of computational memory. This problem can be mitigated by finding the smallest number of Gaussian functions in the approximation that results in an acceptable level of accuracy.

While these proceedings focused on the sum of Gaussian functions approximation, the numerical accuracy of other approximations and numerical integrals was verified during the implementation of the new family of microscopic functionals. However, other sources of uncertainty remain present and need to be quantified and propagated. For example, beyond mean-field effects near closed shell nuclei, which are not described within the single-reference \glsxtrprotectlinksEDF approach used here, can clearly be seen in figure 5 of [12]. This shows a clear and systematic discrepancy between experimental values and the \glsxtrprotectlinksHFB approach and tools for modeling this discrepancy need to be implemented. An approach based on Bayesian Neural Networks has been proposed for this very particular type of discrepancy in simpler mass models like the liquid drop model [20], and a similar approach can be directly applied to density based calculations within the \glsxtrprotectlinksHFB formalism.

5 Acknowledgments

Support for this work was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Nuclear Physics. It was partly performed under the auspices of the US Department of Energy by the Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. Computing support for this work came from the Lawrence Livermore National Laboratory (LLNL) Institutional Computing Grand Challenge program.