Effective theory for low-energy nuclear energy density functionals

J. Dobaczewski, K. Bennaceur, F. Raimondi

July 17, 2012

Abstract:

We introduce a new class of effective interactions to be used within the energy-density-functional approaches. They are based on regularized zero-range interactions and constitute a consistent application of the effective-theory methodology to low-energy phenomena in nuclei. They allow for defining the order of expansion in terms of the order of derivatives acting on the finite-range potential. Numerical calculations show a rapid convergence of the expansion and independence of results of the regularization scale.

A consistent formulation of the low-energy energy-density-functional (EDF) approach in terms of effective theory is long overdue in nuclear physics. Effective theories are now used in different domains of physics with a tremendous success. They quantify and formalize basic physical intuitions about aspects of the description that are essential, and those that can be treated in terms of controlled and correctible approximations. One spectacular example relates to derivations of nuclear forces based on the chiral perturbation theory augmented by the effective-field-theory treatment of corrective terms [1,2]. Analogous effective approach can also be applied in a much simpler context of quantal one-body dynamics [3].

Successful derivations from chiral effective field theory of relativistic [4,5] and nonrelativistic [6,7] nuclear EDFs were recently accomplished. They lead to local EDFs with density-dependent couplings associated with the underlying pion-exchange interaction. These approaches strive to map out the in-medium nucleonic effects at the two-pion-mass (2 $m_\pi\simeq260$MeV/c$^2$) or Fermi-momentum ($k_F\simeq260$MeV/ $\hbar{}c\simeq1.3$fm$^{-1}$) scales. However, the energy scale of low-energy nuclear phenomena in finite nuclei is much lower. For example, to dissociate a nucleon from a nucleus one has to increase its kinetic energy by $\delta{}E_{\mbox{\rm\scriptsize {kin}}}=\hbar^2{}k_F\delta{}k/M\simeq0.25\hbar{}c\,\delta{}k$, that is, a typical separation energy of $\delta{}E_{\mbox{\rm\scriptsize {kin}}}\simeq8$MeV corresponds to the momentum increase of $\delta{}k\simeq32$MeV/$\hbar {}c$. Moreover, the low-energy EDFs are supposed to describe nuclear excitations and shell-effects at even lower energies of 1MeV or $\delta{}k\simeq4$MeV/$\hbar {}c$ and below.

Therefore, what in the QCD-driven chiral dynamics is the small-energy scale, becomes a short-range high-energy scale of nucleon-nucleon force acting on relatively weakly bound nucleons in nuclei. The question of bridging these two different energy scales is the subject of current intense studies in nuclear physics [8]. In finite nuclei, surface effects decrease the infinite-matter particle binding energies by about a factor of two. In the present study, we aim at building a phenomenological low-energy EDFs, with gradient terms providing for the surface terms, corresponding the nuclear leptodermous expansion [9,10].

Ideas of the effective theory have recently been applied in constructing the EDFs based on higher-order derivatives of nuclear one-body densities [11]. When analyzed in terms of the density-matrix expansion [12], the obtained quasilocal EDFs converge fast with the order of derivatives [13], which is a very encouraging feature of the approach. The quasilocal EDFs are also obtained by the Hartree-Fock averaging of higher-order zero-range pseudopotentials [14].

However, one aspect of the approach developed so far is still missing, namely, the quasilocal EDFs and zero-range pseudopotentials do not provide us with any expansion scale, against which the higher-order derivatives could be classified. In the chiral effective field theory, such a scale is given by the cut-off energy in the momentum-space integrals, and allows for formulating a consistent power-counting scheme of the diagrammatic expansion [15]. In low-energy nuclear physics, an analogous cut-off regularization of the zero-range interactions has already been proposed [16,17], and the dimensional regularization was also studied [18], however, a consistent expansion scheme is not yet available.

In the present study, we introduce the expansion scale by employing the regularized higher-order zero-range (pseudo)potentials - exactly as it was illustrated in the simple examples discussed in Ref. [3]. Among a plethora of various ways to regularize delta interactions, we consider the one based on the Gaussian function,

$$\delta(\bm{r}) = \lim_{a\rightarrow\,0} g_a(\bm{r}) = \lim_{... ...\frac{e^{-\frac{\bm{r}^2}{a^2}}}{\left(a\sqrt{\pi}\right)^3} .$$ (1)

Then, the resulting central two-body regularized pseudopotential reads,

$$\hspace*{-0.5cm} V(\bm{r}'_1,\bm{r}'_2;\bm{r}_1,\bm{r}_2) =$$ (2)

$$\sum_{i=1}^4 \hat{P}_i \hat{O}_{i}(\bm{k},\bm{k}')\delta(\bm{r}'_1-\bm{r}_1)\delta(\bm{r}'_2-\bm{r}_2) g_a(\bm{r}_1-\bm{r}_2) ,$$

where $\bm{k}=\frac{1}{2i}(\bm{\nabla}_1-\bm{\nabla}_2)$ and $\bm{k}'=\frac{1}{2i}(\bm{\nabla}'_1-\bm{\nabla}'_2)$ are the standard relative-momentum operators, and the Wigner, Bartlett, Heisenberg, and Majorana terms are given by the standard spin and isospin exchange operators, $\hat{P}_1\equiv 1$, $\hat{P}_2\equiv\hat{P}_\sigma$, $\hat{P}_3\equiv -\hat{P}_\tau$, $\hat{P}_4\equiv -\hat{P}_\sigma\hat{P}_\tau$. The relative-momentum operators are the building blocks of the scalar higher-order derivative operators $\hat{O}_{i}(\bm{k},\bm{k}')$. Generalization of the pseudopotential to spin-orbit and tensor terms are straightforward [14], and in the present study are not discussed.

In our definition of the pseudopotential (2), operators $\hat{O}_{i}(\bm{k},\bm{k}')$ act on the product of the locality delta functions $\delta(\bm{r}'_1-\bm{r}_1)\delta(\bm{r}'_2-\bm{r}_2)$ and regularized delta function $g_a(\bm{r}_1-\bm{r}_2)$. Derivatives of the delta functions have to be understood in the usual sense of distributions, namely, when the pseudopotential (2) is inserted into a matrix element, or when it is averaged with respect to densities, the integration by parts transfers all derivatives onto wave functions or densities. After the integration by parts, the locality deltas reduce the integral to two variables, $\bm{r}_1$ and $\bm{r}_2$, only.

To give a specific example, up to the second-order, that is, up to the next-to-leading-order (NLO) expansion, operators $\hat{O}_{i}(\bm{k},\bm{k}')$ read

$$\hat{O}_{i}(\bm{k},\bm{k}') = T^{(i)}_0 +\frac{1}{2}T^{(i)}_1 \left(\bm{k}^2+\bm{k}'^2\right)+T^{(i)}_2\bm{k}\cdot\bm{k}',$$ (3)

where $T^{(i)}_k$ are the channel-dependent coupling constants. In the limit of $a\rightarrow0$, Eq. (3) links the regularized pseudopotential (2) to the standard zero-range Skyrme interaction [19].

Before considering the general pseudopotential (2), let us first assume that the differential operators $\hat{O}_{i}(\bm{k},\bm{k}')$ depend only on the sum of relative momenta, $\hat{O}_{i}\left(\bm{k}+\bm{k}'\right)$. At NLO, this corresponds to equal second-order coupling constants, $T^{(i)}_1=T^{(i)}_2$. Such particular differential operators commute with the locality deltas $\delta(\bm{r}'_1-\bm{r}_1)\delta(\bm{r}'_2-\bm{r}_2)$, and thus can be applied directly onto the regularized delta $g_a(\bm{r}_1-\bm{r}_2)$. In such a case, pseudopotential (2) reduces to a simple local potential

$$V(\bm{r}) = \sum_{i=1}^4 \hat{P}_i V_{i}(\bm{r}), = \sum_{i=1}^4 \hat{P}_i \hat{O}_{i}(\bm{k}) g_a(\bm{r}) ,$$ (4)

where $V_{i}(\bm{r})=\hat{O}_{i}(\bm{k}) g_a(\bm{r})$ are functions of the relative distance $\bm{r}=\bm{r}_1-\bm{r}_2$. Moreover, since $\hat{O}_{i}(\bm{k})$ are scalar differential operators, the potentials must have forms of power series of Laplacians $\Delta$ in $\bm{r}$, that is,

$$V_{i}(\bm{r}) = \sum_{n=0}^{n_{\mbox{\rm\scriptsize {max}}}} V^{(i)}_{2n} \Delta^n g_a(\bm{r}) ,$$ (5)

where $V^{(i)}_{2n}$ are the coupling constants at order $2n$.

Within the spirit of the effective theory, the series in Eq. (5) is cut at a given $n_{\mbox{\rm\scriptsize {max}}}$. Then, at a given fixed scale $a$, coupling constant $V^{(i)}_{2n}$ have to be adjusted to data. If the expansion is relevant, the results should only weakly depend on the regularization scale $a$, provided it is neither too large, when the expansion scale intrudes into the physical scale of low-energy phenomena, nor too small, where unresolved high-energy features are incorrectly described. At any given order one can pick a specific scale $a$ where the results agree best with the data.

Before such a complete program is carried out, in the present exploratory study we use the fact that a successful parametrization of the local potential already exists in the form of the Gogny interaction [20,21], whose finite-range part can be written as

$$G(\bm{r}) = \sum_{i=1}^4 \hat{P}_i G_i(\bm{r}) = \sum_{i=1}^4 \hat{P}_i \sum_{k=1,2} G^{(i)}_k g_{a_k}(\bm{r}) ,$$ (6)

where Gaussians with two ranges, $a_1=0.7$fm and $a_2=1.2$fm, have been employed. Therefore, instead of adjusting the regularized potential to data, we first attempt to derive it directly from the Gogny interaction. In the present study, we focus on the finite-range part of the Gogny interaction and we keep its zero-range spin-orbit and density-dependent terms unchanged. We then use the results obtained for the original Gogny interaction as metadata, wherewith we compare those obtained for the derived regularized potentials.

At NLO ( $n_{\mbox{\rm\scriptsize {max}}}=1$), the eight coupling constants, $V^{(i)}_0$ and $V^{(i)}_2$, of the local regularized potential (4) can be easily obtained from the eight parameters of the central Gogny force, $G^{(i)}_1$ and $G^{(i)}_2$. Indeed, at this order, potentials (5) are combinations of a Gaussian and Gaussian multiplied by $\bm{r}^2$, which span similar sets of functions as those of two Gaussians with two different ranges.

Below we use a prescription that can be applied at any order $n_{\mbox{\rm\scriptsize {max}}}$, and which is based on global characteristics of the potentials. Namely, we determine coupling constants $V^{(i)}_{2n}$ in Eq. (5) by requiring that the lowest moments of both potentials are equal, that is,

$$M_{2m}^{(i)} \equiv \int \bm{r}^{2m} G_i(\bm{r}) {\rm d}^3 \bm{r} = \int \bm{r}^{2m} V_i(\bm{r}) {\rm d}^3 \bm{r} ,$$ (7)

for $m=0,1,\ldots,n_{\mbox{\rm\scriptsize {max}}}$. This conditions gives the coupling constants of the regularized potential in simple analytical forms,

$$V^{(i)}_{2n} = \sum_{m=0}^n \left(-\frac{a^2}{4}\right)^{n-m} \frac{M_{2m}^{(i)}}{(n-m)!(2m+1)!}$$

$$= \frac{1}{4^n n!}\sum_{k=1,2} G^{(i)}_k\left(a^2_k-a^2\right)^n ,$$ (8)

where $G^{(i)}_k$ and $a_k$ are the parameters of the Gogny interaction, as in Eq. (6).

To determine properties of finite nuclei, in the present study we used the EDF approach for regularized potential defined in Eq. (4), along with the coupling constants (8) derived for the D1S Gogny force [21]. Numerical calculations were performed by implementing regularized potentials in the code HFODD (v2.53r) [22,23], which was adapted to solving self-consistent equations for finite-range interactions.

In Fig. 1, we present results of calculations performed for eight doubly magic nuclei: $^{16}$O, $^{40}$Ca, $^{48}$Ca, $^{56}$Ni, $^{78}$Ni, $^{100}$Sn, $^{132}$Sn, and $^{208}$Pb. In each case, we show relative deviations of binding energies (left panels) and rms radii (right panels) as compared to values obtained with the original Gogny interaction. One can see that with increasing order of expansion, the convergence is very rapid; indeed, for every next order, the obtained relative deviations decrease by about a factor of four, and at N$^3$LO they are well below 1%.

Figure 1: (Color online) Deviations (in percent) of the EDF binding energies [left panels (a), (c), and (e)] and rms radii [right panels (b), (d), and (f)], calculated for the regularized potential (4), relative to values obtained with the original Gogny interaction. Top, middle, and bottom panels correspond to the NLO, N$^2$LO, and N$^3$LO expansions, that is, to $2n_{\mbox{\rm\scriptsize{max}}}$=2, 4, and 6, respectively. The legend gives the atomic masses of doubly magic nuclei listed in the text. At NLO, missing points correspond to large values of $a$, where converged self-consistent solutions could not be obtained.

We note here that the results presented in Fig. 1 were obtained without any adjustment of parameters. However, one can expect that such small adjustments can easily correct for the smooth trends observed in Fig. 1. This is so because binding energies of nuclei, which result form cancellations of the kinetic and potential energies, are extremely sensitive to parameters of the interaction. For example, we checked that modifications of the Gogny density-dependent term of the order of 0.05% are able to bring the N$^3$LO mass of $^{208}$Pb exactly to the Gogny value.

The regularized potentials proposed in this study (with or without the restricted dependence on the relative momenta introduced above), will later be used in extensive optimization procedures - similarly as it has been done for other standard functionals used within the EDF approaches. Nevertheless, even without such extensive adjustments, basic features of the obtained effective theory can already be illustrated by studying scaling properties of the results. To this end, in Fig. 2 we show deviations of binding energies and radii relative to those obtained for $^{208}$Pb. One can see that the results scaled in this way are almost independent of the regularization scale, which is an extremely encouraging result, consistent with the effective-theory principles.

Figure 2: (Color online) Same as in Fig. 1, but for the deviations relative to those obtained for $^{208}$Pb. For clarity, near $a\simeq 1$fm, where the $^{208}$Pb normalization values cross zero (cf. Fig. 1), several points were removed from the plot.

Effective theories are built to provide us with expansions that converge better when the bound-state energy or scattering energy is lower. Spectacular examples of such a convergence pattern exist for weakly bound electron states in quantum mechanics of a one-body Coulomb problem [3] or for phase shifts of the two-body nucleon-nucleon scattering [2]. Convergence of an effective theory can be best visualized by the so-called Lepage plot [3], which shows the dependence of the error in the description of a given observable on energy. Such error should decrease with decreasing energy and with increasing order of expansion.

In many physical systems, scale of energy is inverted with respect to the scale of distance, that is, small distances correspond to large energies and vice versa. In finite nuclei, the situation is more complex, because the average energy per particle and average internucleon distance are almost constant throughout the mass chart, which is a simple consequence of the saturation of nuclear forces and near-constancy of density inside nuclei. Should then the effective theories work better in light (small) nuclei or in heavy (large) ones?

Good arguments in favor of the second option are provided by the density-matrix expansion (DME) [12], which uses the fact that the range of nuclear forces is smaller than characteristic distances at which densities in nuclei vary, see Refs. [24,13] and references cited therein, and by the ideas of the leptodermous expansion [9,10]. Indeed, since the largest variations of densities occur at nuclear surfaces, effective theories should work best in large nuclei, where the surface region is relatively small with respect to the nuclear bulk.

To study this problem, in Fig. 3 we show reduced deviations of binding energies, $\Delta b=\Delta B/A$, and radii, $\Delta r=\Delta R/A^{1/3}$, in function of the particle number $A$. This allows us to compare results obtained for different nuclei in the same scale. Since the relative deviations are almost independent of the regularization scale $a$ (cf. Fig. 2), here we plot them at $a=0.85$fm.

Figure 3: (Color online) Same as in Fig. 2, but for the deviations of binding energies scaled by the numbers of particles, $\Delta b=\Delta B/A$, and those of radii scaled by the nuclear sizes, $\Delta r=\Delta R/A^{1/3}$. Note that the N$^2$LO results almost perfectly hide behind those obtained at N$^3$LO.

We see that the results obtained for nuclei beyond about $A=48$ scale with mass as those for $^{208}$Pb. From the point of view of the effective theory and proposed expansion, one can interpret this result as a high similarity of these nuclei, that is, their sizes are not different enough to induce significant differences in their convergence properties, and they all would appear at the same point of the Lepage plot. We also see that properties of lighter nuclei, with $A\leq 48$, converge slightly better than those of the heavier ones. However, only a future complete analysis, performed for the entire mass table and for optimized functionals, can here provide conclusive results.

We complete our analysis by discussing values of coupling constants (8) in the so-called natural units, see Refs. [25,26] and references cited therein. Since the EDF generated by our regularized potential contains only terms quadratic in densities, in the natural units the coupling constants have the following dimensionless values,

$$v^{(i)}_{2n}= f^2 \Lambda^{2n} V^{(i)}_{2n},$$ (9)

where $f^2$ is an overall normalization factor and $\Lambda^2$ is a characteristic scale of the derivative operator $\Delta$.

Coupling constants $V^{(i)}_{2n}$ (8) are plotted in Fig. 4. Their almost-linear $n$-dependence, which we see in the logarithmic scale, allows us to derive the slope of this dependence that corresponds to the value of $\Lambda \simeq 700$MeV/ $\hbar{}c\simeq3.5$fm$^{-1}$ (dashed line). For our regularized potentials at N$^3$LO, we thus obtain a similar scale $\Lambda$ as that characterizing the zero-range interactions at NLO ( $2n_{\mbox{\rm\scriptsize {max}}}$=2) [25,26]. Next, we fix the normalization factor $f$ so as to make the values of dimensionless coupling constants $v^{(i)}_{2n}$ of the order of unity, which gives $f\simeq35$MeV/ $(\hbar{}c)^{3/2}\simeq1$/(77MeV$^{1/2}$fm$^{3/2}$). As we see, this value is significantly smaller than the pion decay constant $f_\pi\simeq93$MeV/ $(\hbar{}c)^{3/2}$, which normalizes the EDFs derived directly from the chiral effective field theory [25].

Figure 4: (Color online) Coupling constants $V^{(i)}_{2n}$ derived at $a=0.85$fm from the Gogny interaction, as dictated by Eq. (8). The dashed line illustrates the slope of the $n$-dependence of the coupling constants, plotted here for $\Lambda \simeq 700$MeV/$\hbar {}c$.

Finally, for the derived values of $\Lambda$ and $f$, in Fig. 5 we show the dependence of the dimensionless coupling constants $v^{(i)}_{2n}$ on the regularization scale $a$. We see that for $n=1$, 2, and 3 they have values of the order of unity at all scales, although the fourth-order coupling constants ($n=2$) cross zero around $a=0.95$fm. The zero-order coupling constants ($n=0$) appear to be less natural, with values of the order of 0.1 (cf. also Fig. 4). However, we should keep in mind that in the present study we derived the coupling constants from the Gogny interaction, which has a predefined scale of the order of 1fm. Therefore, we can expect that the future adjustments of the coupling constants directly to data may lead to even weaker scale-dependence.

Figure 5: (Color online) Coupling constants $v^{(i)}_{2n}$ in natural units (9) as functions of the regularization scale $a$, plotted for $n=0$-3 in panels (a)-(d).

In summary, in this study we introduced a new class of energy density functionals that are based on pseudopotentials related to regularized delta interactions. We used the regularization scheme employing the standard Gaussian form, which introduces a regularization scale $a$ corresponding to the range of the interaction. Different orders of expansion then correspond to different orders of derivatives used in the pseudopotential.

In our opinion, future prospects for using the proposed regularized (pseudo)potentials are high. First, similarly as in the simple one-body examples [3], they may present better convergence properties than similar expansions based on the zero-range interactions. Second, they allow for convergent summations of contributions from high single-particle momenta, which is not the case for zero-range interactions. And third, they allow for formulating a consistent expansion in terms of the orders of derivatives, with the convergence properties gauged against the regularization scale. This last feature is unique among all the EDF approaches based on zero-range and finite-range interactions developed so far, and gives us a potential of building an order-by-order correctible theory.

Interesting comments by Dick Furnstahl and Witek Nazarewicz are gratefully acknowledged. This work has been supported in part by the Academy of Finland and University of Jyväskylä within the FIDIPRO programme and by LEA COPIGAL Polish-French Collaboration Agreement. We acknowledge the CSC - IT Center for Science Ltd, Finland for the allocation of computational resources.