Convergence of density-matrix expansions for nuclear interactions

B.G. Carlsson and J. Dobaczewski

July 28, 2010

Abstract:

We extend density-matrix expansions in nuclei to higher orders in derivatives of densities and test their convergence properties. The expansions allow for converting the interaction energies characteristic to finite- and short-range nuclear effective forces into quasi-local density functionals. We also propose a new type of expansion that has excellent convergence properties when benchmarked against the binding energies obtained for the Gogny interaction.

The description of finite many-body systems is a continuous exciting challenge in physics. In view of the exploding dimensionality of many-particle Hilbert spaces, direct approaches are however often untractable. Instead, density functional theory (DFT) [1] provides an alternative by showing that a much simpler description in terms of the density alone is in principle possible. The main challenge is finding or deriving the best density functionals that are able to consistently describe large classes of physical systems. In this work we address a related question, namely to what extent the physics of the exact nonlocal exchange potential can be captured by an approximate density functional? This question is not only of fundamental importance but also of great practical interest, and attempts to solve it date back to the early 1950s, involving the Slater as well as the optimized effective potential approximations [1].

For nuclear-physics applications, Negele and Vautherin developed an alternative method known as the density-matrix expansion (DME) [2], which takes advantage of the short-range of the nuclear two- and three-body forces. It is particularly convenient as it was originally constructed to generate quasi-local density functionals [3] of a form similar to the successful phenomenological Skyrme functionals used extensively today [4]. The method have recently gone through a revival [5,6,7,8,9] and appears to be a promising starting point for the construction of ab initio nuclear DFT models [7,8]. The idea is to start from effective forces, derived from first principles, and then transform the expression for the interaction energies into quasi-local density functionals using the DME. Because the starting point is an effective force, this method is able to take correlations into account in a way that goes much beyond the Hartree-Fock treatment of the bare interaction.

In conjunction with these ab initio methods, there are empirical efforts to improve nuclear density functionals, such as the systematic method we recently proposed to construct nuclear density functionals using higher-order derivatives of densities [10]. In such an approach, the DME is helpful as it can serve as a method of classifying different terms according to their importance, similarly to what is available in the effective-field-theory power-counting scheme [11].

In this Letter, we explore the ideas of the DME to transform the direct and exchange energies resulting from the finite-range Gogny interaction into quasi-local density functionals. We then extend these methods beyond the second order expansions used so far to make the first tests of convergence and accuracy in function of the expansion order.

We begin our analysis by recalling that the nuclear one-body density matrix $\rho\left(\vec{r}_{1}\sigma_{1}\tau_{1},\vec{r}_{2}\sigma_{2}\tau_{2}\right)$, which depends on the space-spin-isospin nucleon coordinates $\vec{r}\sigma\tau$, can be separated into four nonlocal spin-isospin densities $\rho_{v}^{t}\!\left(\vec{r}_1,\vec{r}_2\right)$ as

$$\rho\left(\vec{r}_{1}\sigma_{1}\tau_{1},\vec{r}_{2}\sigma_{2}\tau_{2}\right)=$$

$$~~~{\textstyle{\frac{1}{4}}}\sum_{vt}\left(\sqrt{3}\right)^{v+t} ... ...au_{2}}, \rho_{v}^{t}\!\left(\vec{r}_1,\vec{r}_2\right)\right]^{0}\right]_{0} ,$$ (1)

where we sum over the spin-rank $v=0$ and 1 (scalar and vector) and isospin-rank $t=0$ and 1 (isoscalar and isovector) spherical tensors of the Pauli matrices, $\sigma_v$ [12] and $\tau^t$, respectively. We denote vector (isovector) coupling by the standard square brackets with subscripts (superscripts) giving the values of the total spin (isospin).

For a local interaction, the energy corresponding to the direct (Hartree) term depends on the local densities $\rho_{v}^{t}\!\left(\vec{r}\right)\equiv\rho_{v}^{t}\!\left(\vec{r},\vec{r}\right)$ in Eq. (1) as

$${\cal E}^{\mbox{\rm\scriptsize {int}}}_{\mbox{\rm\scriptsize {dir}}} =\!\!\! {\textstyle{\frac{1}{2}}}\!\!\int\!\! {\rm d}\vec{r}_{1}{\rm d}\v... ...c{r}_{1}\right),\rho_{v}^{t}\!\!\left(\vec{r}_{2}\right)\right]^{0}\right]_{0},$$ (2)

where $V_{vt}^{\mbox{\rm\scriptsize {dir}}}\left(\vec{r}_{1},\vec{r}_{2}\right)$ denote the local direct two-body potentials. As discussed in Ref. [9], for short-range interactions and sufficiently smooth local densities one can employ simple Taylor expansions of densities to obtain the energy density in the form of a gradient expansion. As in Ref. [10], we use the spherical tensor representation of derivatives, which allows for the most economical classification of terms in the expansions.

Introducing the standard average $\vec{R}=\frac{1}{2}\left(\vec{r}_1+\vec{r}_2\right)$ and relative $\vec{r}=\left(\vec{r}_1-\vec{r}_2\right)$ coordinates, we have the following expansion of the local densities,

$$\rho_{v}^{t}\!\left(\vec{R}+\frac{\vec{r}}{2}\right) =\! e^{{\textstyle{\frac{1}{2}}}\vec{r}\cdot\hat{\vec{\nabla}}} \!\!\... ...}}}\vec{r}\cdot\hat{\vec{\nabla}}\right)^{n} \rho_{v}^{t}\!\left(\vec{R}\right)$$

$$=\! \sum_{nL}\,a_{nL}\,r^n\,\left[Y_{L}\left( \hat{r} \right),\hat{D}_{nL}\right]_{0} \rho_{v}^{t}\!\left(\vec{R}\right),$$ (3)

where, for clarity, the spin and isospin projection quantum numbers, $m_v$ and $m_t$, are not shown, $\hat{\vec{\nabla}}$ denotes the gradient operator with respect to variable $\vec{R}$, $\hat{D}_{nL}$ denote the spherical-tensor derivatives of order $n$ and rank $L$ [10], $Y_{L}\left( \hat{r} \right)$ are the standard spherical harmonics, and $a_{nL}$ are simple numerical constants resulting from the recoupling of the Cartesian $n$th order derivatives $\left({\textstyle{\frac{1}{2}}}\vec{r}\cdot\hat{\vec{\nabla}}\right)^{n}$ into the spherical-tensor forms [13].

Table 1: The direct end exchange energies in $^{208}$Pb calculated with the finite-range Gogny D1S interaction [15] in the scalar ($v=0$), vector ($v=1$), isoscalar ($t=0$), and isovector ($t=1$) channels. All energies are in MeV.

		$t=0$		$t=1$
		$v=0$		$v=1$		$v=0$		$v=1$
${\cal E}^{\mbox{\rm\scriptsize {int}}}_{\mbox{\rm\scriptsize {dir}}}$		$-$12294.769		0		373.314		0
${\cal E}^{\mbox{\rm\scriptsize {int}}}_{\mbox{\rm\scriptsize {exc}}}$		$-$594.992		0.629		$-$32.337		0.215

Expanded in this way, and after a partial integration and recoupling needed to make all derivative operators act on the same density, the direct interaction energy of Eq. (2) can be written as a gradient expansion,

$${\cal E}^{\mbox{\rm\scriptsize {int}}}_{\mbox{\rm\scriptsize {dir}}} \simeq\!\! \sum_{mIvt}\int\!\!{\rm d}\vec{R}\, C_{mIv}^{t} \left[ \left[\rho... ...{R}),\left[\hat{D}_{mI},\rho_{v}^{t}(\vec{R})\right]_{v}\right]^{0}\right]_{0}.$$ (4)

In this final form, the coupling constants $C_{mIv}^{t}$ are numbers obtained as moments of the direct potentials up to the order of the expansion $N$. This order is defined by keeping all terms in the product that fulfill the condition $n+n' \leq N$ , where $n$ and $n'$ denote the orders in the expansions of the two densities [compare with Eq. (3)].

To test the quality and convergence properties of this expansion, we generated a set of realistic test densities for doubly magic nuclei by solving the self-consistent equations for the nuclear SLy4 functional [14]. For these test densities, the sums in Eq. (4) were evaluated to different orders by calculating the higher-order derivatives of densities needed using the program HOSPHE (v1.00) [12]. The coupling constants in Eq. (4) were calculated from moments of the finite-range part of the Gogny D1S interaction [15]. Then for the same set of test densities, the exact direct and exchange energies of the Gogny interaction were determined using the capability of treating Gaussian interactions with the code HFODD (v2.45h) [16].

As an example, in Table 1 we show the exact energies calculated in $^{208}$Pb for the four spin-isospin channels. In the ground states of even-even nuclei, the vector terms of the direct energies vanish due to the conserved time-reversal symmetry, and thus our calculations here only test the scalar terms of the expansion in Eq. (4). In Fig. 1 we show percentage deviations between the gradient approximations and the exact results. It is extremely gratifying to see that in light nuclei the sixth-order expansion allows for reaching the precision of 0.1-0.01% or 0.1MeV. Even more important, in each higher order the precision increases by a large factor, which is characteristic to a rapid power-law convergence. One also sees that in heavier, isospin-unsaturated nuclei the precision slightly deteriorates, probably owing to a slower convergence for the isovector densities, which fluctuate more than the isoscalar ones.

Figure 1: (Color online) Precision of the gradient expansions (4) of the direct interaction energies calculated for the finite-range Gogny D1S interaction [15]. The nine nuclei used for the test are $^4$He, $^{16}$O, $^{40,48}$Ca, $^{56,78}$Ni, $^{100,132}$Sn, and $^{208}$Pb. Positive (fourth order) and negative (second and sixth orders) errors are shown with open and closed symbols, respectively.

One can, of course, quite easily calculate the direct terms exactly, but still it is interesting to see that a few of the lowest-order coupling constants appearing in the expansion are able to capture the essential physics contained in the nuclear finite-range force. For the exchange (Fock) terms whereto we turn now, the situation is significantly different, because here the expansion into quasi-local densities dramatically speeds up the calculations.

In order to deal with the exchange term, we propose a new approximation, which we call a damped Taylor (DT) expansion, and we define it in the following way:

$$\rho_{v}^{t}\!\left(\vec{R},\vec{r}\right)\!=\! \bar{\rho}_{... ...\hat{r}\right),\rho_{nLv}^{t}\left( \vec{R}\right)\right]_{0}.$$ (5)

The Taylor expansion of the term in the square brackets is obtained by using $e^{i\vec{r}\cdot\hat{\vec{k}}}=\sum_{m}\frac{1}{m!}\left(i\vec{r}\cdot\hat{\vec{k}}\right)^{m},$ where $\hat{\vec{k}}={\textstyle{\frac{1}{2i}}}\left(\vec{\nabla}_1-\vec{\nabla}_2\right)$ is the standard relative-momentum operator. The quasi-local densities, are defined as

$$\rho_{nLv}^{t}\left( \vec{R}\right) =\left.\hat{K}_{nL} \rho... ...c{r}_2\right) \right\vert _{\vec{r}_{1}=\vec{r}_{2}=\vec{R}} ,$$ (6)

where the higher-order spherical-tensor derivative operators $\hat{K}_{nL}$ [10] are built from the first-order operators $\hat{\vec{k}}$. Because of the simple form of the expansion, explicit formulas for functions $\pi_{nL}^m\left(r\right)$ can be derived and will be presented elsewhere [13].

Figure: (Color online) Comparison of a Taylor expansion and the DT expansion (with $\bar{\rho}_{v}^{t}=0$ and $a=4/k_F$), applied to the SNM nonlocal density in second, fourth, and sixth orders. The figure is drawn for $k_F^0=1.35$ fm$^{-1}$ corresponding to the equilibrium density of SNM [15,9]. The shaded regions indicate the integration ranges needed to obtain 1, 0.1, 0.01 and 0.001 % accuracy for the SNM exchange energy with the Gogny D1S interaction.

The DT expansion has the potential to be more accurate than the original Negele-Vautherin (NV) expansion [2] because it avoids the approximation of angular averaging the density matrix in the nonlocal $\vec{r}$ direction. It is generated around a model density $\bar{\rho}_{v}^{t}$ and is also damped by a Gaussian factor with width $a$. In the limit of vanishing model density and infinite damping distance, the DT becomes an ordinary Taylor expansion. This Taylor expansion applied to the uniform system of symmetric nuclear matter (SNM) is illustrated in Fig. 2. It works well in a region around $\vec{r}=0$ but diverges for larger relative distances. These bad asymptotics can however be cured by using a finite damping width $a$, as also demonstrated in Fig. 2. By employing the condition that it should work equally well in SNM for all values of Fermi momentum $k_F$, we estimate this width as $a=4/k_F$. Assuming that the densities of real nuclei have similar behavior, we retain this estimate of $a$ in the calculations presented below.

The model density can be chosen in different ways by the requirement that the truncated expansions must become exact in certain limits. The model used in the following leads to exact values in SNM, and is obtained by assuming that the density in each point $\vec{R}$ has the SNM nonlocal dependence, $\bar{\rho}_{v}^{t}\!\left(\vec{r}_1,\vec{r}_2\right) =\rho_{v}^{t}\!\left(\vec{R}\right)\frac{3j_{1}\left(k_{F}r\right)}{k_{F}r}.$ For large relative distances, when the expansion around $r=0$ can no longer be trusted, the bracketed term in Eq. (5) is damped to zero and the asymptotic behavior is then given entirely by this model part. In this way the parameter $a$ interpolates between placing trust in the model part or in the expansion part, and the optimal balance is determined by the range in which the expansion converges, as estimated above.

By inserting the DT expansion of Eq. (5) into the expression for the exchange energy, we obtain the quasi-local expansion

$$\!\!{\cal E}^{\mbox{\rm\scriptsize {int}}}_{\mbox{\rm\script... ... \rho_{nLvJ}^{t}\!\left(\vec{R}\right)\right]^{0}\right]_{0} ,$$ (7)

where the quasi-local densities $\rho_{nLvJ}^{t}\!\left(\vec{R}\right)$, which were introduced in Ref. [10], are equal to those of Eq. (6) with the angular momenta $L$ and $v$ coupled to the total angular momentum $J$. The coupling-constants $C_{nLvJ}^{t,n'}$ are given by the integrals of the local exchange two-body potentials $V_{vt}^{\mbox{\rm\scriptsize {exc}}}\left(\vec{r}_{1},\vec{r}_{2}\right) =V_{vt}^{\mbox{\rm\scriptsize {exc}}}\left(r\right)$ and functions $\pi_{nL}^m\left(r\right)$ of Eq. (5). They depend on $k_F$ through the model density and damping width, which in the standard local density approximation leads to a density dependence $k_F=(3\rho_0^0/2\pi^2)^{1/3}$ [2].

Figure 3: (Color online) The RMS deviations between the exact and approximate exchange energies of Eq. (7), calculated for the nine nuclei listed in the caption of Fig. 1. The four panels show results obtained in the four spin-isospin channels labeled by $V_{vt}$. Results obtained in the present work (DT) are compared with those corresponding to the NV [2] and PSA [7] expansions.

In Fig. 3 we show results obtained with the DT expansion compared to the exact values, similar to those listed in Table 1. The comparison is made by considering the RMS deviations between the exact and approximate exchange energies, and is also shown for the NV [2] and phase-space averaging (PSA) [7] expansions. As seen in this figure, the NV Bessel-function expansion used in the scalar channels works very well, whereas the alternative method used for the vector channels [see Eq. (4.23) in Ref. [2]] does not improve when evaluated beyond the second order.

We employ the recently proposed PSA method in the INM-DME version [7] and generalize it to higher orders by incorporating the model density in the same way as in Eq. (5). Other modifications are probably also possible and could be studied as well. At second order, our prescription reproduces that proposed in Ref. [7] with the PSA performed over SNM. Fig. 3 shows that this is the best expansion at second order, where the total RMS deviation for the exchange energy is 14.09 MeV. When used for the vector part, the PSA expansion is similar to the NV expansion and also has problems when evaluated beyond second order.

As seen in Fig. 3, the DT expansion exhibits excellent convergence properties both for the vector and scalar parts, with every next order improving the results by large factors. Note that in the vector channels, the exact exchange energies of Table 1 are quite small, so at second order their relative description is poor. In contrast, at higher orders their description becomes rather accurate. The results obtained for the DT expansion are not so sensitive to the value of the damping width $a$, and values in the range of $a\approx3/k_F$ - $5/k_F$ all lead to improvements over the NV expansion.

The total RMS deviations obtained in the DT case for the exchange energies at fourth and sixth orders are 1.54 and 0.36MeV, respectively [17]. Calculations using a fixed value of $k_F^0=1.35$fm$^{-1}$ give results that are not significantly different [13]. They lead to substantial simplifications, because then the density dependencies of coupling constants can be neglected.

In summary, we presented the first analysis of density-matrix expansions extended beyond second order, with derivatives of densities included up to sixth order. We also proposed a new type of expansion that has excellent convergence properties. The results demonstrate that its possible to construct expansions into quasi-local density functionals that converge for the description of low-energy nuclear observables. This also facilitates the difficult adjustments of phenomenological higher-order nuclear functionals by providing inital guesses for free parameters as well as a method to derive them.

We thank Johan Grönqvist for useful discussions. This work is supported in part by the Academy of Finland and the University of Jyväskylä within the FIDIPRO programme, and by the Polish Ministry of Science under Contract No. N N202328234.