Fitting Skyrme functionals using linear response theory

The formalism of the linear response (LR) theory has been presented in a recent series of papers written by the Lyon group [3–5] for the case of a general Skyrme functional [6] including spin–orbit and tensor interactions. Such a formalism can be adopted as a useful tool for detecting the instabilities that plague Skyrme functionals [1, 2, 7].

In this paper, we present an improved fitting protocol to build Skyrme functionals that, for the first time, includes the stability of symmetric nuclear matter (SNM) up to a given density as an additional constraint. Such a method is more powerful than the standard constraints on the Landau parameters, since the latter represent the long-wavelength limit and are therefore not able to detect finite-size instabilities.

An example of this kind of instability has been revealed by Lesinski et al [7] and was shown to be related to large values of C^Δρ₁, the coupling constant of the ρ₁Δρ₁ term in the Skyrme functional [6]. This coupling constant, if not properly constrained during the fitting procedure, can drive the merit function χ² used to adjust the parameters of the interaction in a region of matter fluctuations that give unphysical properties for the ground state of some nuclei. Since this gradient term does not contribute to the Landau parameters, this kind of instability is only detectable by the LR method.

This paper is organized as follows. In section 2, we recall the main properties of a Skyrme functional and in section 3 we briefly recall the LR formalism; in section 4, we present the results of our new fitting method and, finally, we present our conclusions in section 5.

The density-dependent Skyrme interaction considered in this paper has the standard form

$$\begin{eqnarray} &&v(\mathbf{R},\mathbf{r}) = t_{0}(1+x_{0}\hat{P}_{\sigma})\delta(\mathbf{r})+ \textstyle\frac{1}{6}\,t_{3}(1+x_{3}\hat{P}_{\sigma}) \rho_{0}^{\alpha}(\mathbf{R})\delta(\mathbf{r})\nonumber\\ &&\quad+\, \textstyle\frac{1}{2}t_{1}(1+x_{1}\hat{P}_{\sigma}) [ {\mathbf{k}\smash{'}}^2\delta(\mathbf{r}) +\delta(\mathbf{r})\mathbf{k}^{2}] + t_{2}(1+x_{2}\hat{P}_{\sigma}) \mathbf{k}'\cdot \delta(\mathbf{r})\mathbf{k}\nonumber\\ &&\quad +\, \mathrm{i}\,W_{0}(\hat{\boldsymbol{\sigma}}_{1} +\hat{\boldsymbol\sigma}_{2}) \cdot [ \mathbf{k}'\times \delta(\mathbf{r})\mathbf{k}], \end{eqnarray} \tag{ 1 }$$

where we used r ≡ r₁ − r₂ and $\mathbf {R}=\frac {1}{2}\left (\mathbf {r}_{1}+\mathbf {r}_{2}\right )$ for the relative and center-of-mass coordinates, respectively. We also defined $\hat {P}_{\sigma }$ as the spin exchange operator, while $\mathbf {k}=-\frac {\mathrm i}{2}(\nabla _{1}-\nabla _{2})$ is the relative momentum acting on the right and k' its complex conjugate acting on the left. Finally, ρ₀(R) is the scalar–isoscalar density of the system. The energy density functional (EDF) is derived by keeping all the terms in the particle–hole channel except the Coulomb exchange one that is treated using the Slater approximation. The center-of-mass correction is approximated by its one-body part.

The random phase approximation (RPA) response function in each channel (α) = (S,M,I) ≡  (spin, projection of the spin and isospin) is obtained by solving the Bethe–Salpeter equations for the correlated Green functions G^(α)_RPA written as

$$\begin{eqnarray} G_{\mathrm{RPA}}^{\alpha}(q,\omega,\mathbf{k}_1)&=&G_{\mathrm{HF}}(q,\omega,\mathbf{k}_1)+G_{\mathrm{HF}}(q,\omega,\mathbf{k}_1)\nonumber\\ & &\times\sum_{\alpha'}\int\!\frac{\mathrm d^3 k_2}{(2\pi)^3}\, V_{\mathrm{ph}}^{\alpha;\alpha'}(q,\mathbf{k}_1,\mathbf{k}_2)\nonumber\\ &&\times G^{\alpha'}_{\mathrm{RPA}}(q,\omega,\mathbf{k}_2), \end{eqnarray} \tag{ 2 }$$

where G_HF(q,ω,k₁) is the unperturbed particle–hole propagator, and V^α;α'_ph(q,k₁,k₂) is the residual interaction obtained performing the second functional derivative. See [3, 4, 8] for more details. Once equation (2) is solved, we can calculate the response function of the system with

$$\begin{equation} \chi^{\alpha}(\omega,q) =4\int\!\frac{\mathrm d^3k}{(2\pi)^{3}}\, G_{\mathrm{RPA}}^{\alpha}(q,\omega,\mathbf{k}), \end{equation} \tag{ 3 }$$

where ω is the transferred energy and q is the transferred momentum. As usual we chose the momentum q oriented along the z-axis of the reference frame. We also use natural units ℏ = c = 1 to simplify the notations. When an instability appears in the infinite system, it manifests itself by a pole in the response function χ^α(ω,q) at zero energy.

In figure 1 (left panel), we show the positions of the poles in SNM for the Skyrme SkP [9] functional. In figure 1 (right panel), we show the imaginary part of the response function in the channel (S,M,I) = (0,0,1) for the same functional and for ρ_pole = 0.205 fm⁻³ and q = 2.88 fm⁻¹.

Zoom In Zoom Out Reset image size

Lesinski et al [7] showed that such a functional is not stable: performing high-precision Hartree–Fock calculations (HF) using a sufficient number of iterations, we observe a strong oscillating isovector density ρ₁; such oscillations increase at each iteration and prevent the code from converging. The same conclusion was drawn for other Skyrme functionals, for example the LNS [10] one. Such unphysical configurations were discovered to be related to an excessive contribution of the term proportional to ρ₁Δρ₁ in the functional.

Performing extensive HF calculations for even–even [2] and odd–even nuclei [11], we found that there exists a relation between the position of the poles of the response function in SNM and the instabilities in finite nuclei. To guarantee that a nucleus will not manifest finite-size instabilities, we have imposed the empirical constraint

$$\begin{equation} \rho_{\mathrm{pole}}>1.4\,\rho_{\mathrm{sat}}, \end{equation} \tag{ 4 }$$

where ρ_sat is the saturation density of SNM. Such a relation has to be satisfied for all the spin/isospin channels. The only instability we allow to exist in the density region ρ < 1.4 ρ_sat is the spinodal instability, which manifests in the (0,0,0) channel at low densities and low transferred momenta and characterizes the liquid–gas phase transition of nuclear matter. Using such a criterion allows us to quickly test large sets of Skyrme functionals. This is a great advantage over self-consistent cranking calculations [12] and systematic quasi-particle blocking calculations [13] used for revealing the instabilities in spin channels.

Usually, the time-odd terms of the functional are not directly constrained by observables. Indirect constraints may come from the fact that some functionals have been constructed with the requirement that the time-even and -odd terms lead to Landau parameters satisfying, respectively, the conditions

$$\begin{equation} \frac{F^{}_{\ell}}{2\ell+1}>-1\,,~~~~\frac{G^{}_{\ell}}{2\ell+1}>-1 \end{equation} \tag{ 5 }$$

for ℓ = 0, 1 and equivalent ones for F'_ℓ and G'_ℓ. These conditions are used to prevent the appearance of global spin and/or isospin instabilities, but do not prevent the formation of finite-size instabilities due to gradient terms such as ρ₁∇ρ₁, s_tΔs_t and (∇·s_t)² for isospin t = 0,1.

An example of such a spin instability can been found in the SLy5 [14] functional. In figure 2, we show the position of the poles of the response function in SNM. We note that there is a pole in the channel (1,1,0) close to saturation density that does not respect the criterion given in equation (4). Performing time-dependent-Hartree–Fock calculations, Fracasso et al [15] have shown that the functional SLy5 is indeed unstable against spin-dipole excitations. The same instability also manifests when performing cranking calculations of finite nuclei [12].

Zoom In Zoom Out Reset image size

Using the results discussed in the previous section, we now present a new fitting protocol to build stable Skyrme functionals. We use a modern version of the Saclay–Lyon fitting protocol similar to the one presented in the recent paper of Washiyama et al [16]. On top of such a protocol, we add the constraint given by equation (4). This constraint is enforced at each iteration of the fit and does not significantly change the execution time since, as already discussed in [4], the expressions for the response functions are analytical.

Using the SLy5 functional as a starting point, we decided to build a new functional with similar properties, but without instabilities in the spin channel. The resulting functional has been called SLy5* and its coupling constants are reported in table 1. In the same table, we compare nuclear matter properties at the saturation point for the new SLy5* and the original SLy5, showing that they are all essentially equal but for the Thomas–Reiche–Kuhn enhancement factor κ_v which is larger with SLy5*. The value obtained with SLy5* is in agreement with the empirical value [20] close to κ^emp_v ≈ 0.4.

4.1. Equations of state

In order to compare the isoscalar and isovector properties of the functionals SLy5 and SLy5*, we have calculated the equation of state of SNM and pure neutron matter (PNM) as well as the neutron effective mass in these two extreme media.

In the upper panel of figure 3, we show the equations of state for SNM and PNM. We observe very good agreement between the results given by the two functionals. In SNM the two curves stay on top of each other and in PNM this is the case for up to two to three times the saturation density.

Zoom In Zoom Out Reset image size

In the lower panel of figure 3, we report the evolution of the neutron effective mass in SNM and PNM. In SNM, the curves that give the neutron effective mass m*/m as a function of the density stay on the top of each and cannot be distinguished. In the PNM case, the neutron effective mass obtained from SLy5 is systematically lower than the effective mass in SNM.

One can also observe that SLy5* almost gives the same evolution for the neutron effective mass in SNM and in PNM. The origin of this property can be clarified by writing the expressions for the neutron effective mass in the two media in terms of coupling constants [7]

$$\begin{equation} \frac{\hbar^{2}}{2m^{*}_{n}}=\frac{\hbar^{2}}{2m}+C^{\tau}_{0}\rho\phantom{++C^{\tau}_{1}} \;\;\;\;\; \textrm{(SNM),} \end{equation} \tag{ 6 }$$

$$\begin{equation} \frac{\hbar^{2}}{2m^{*}_{n}}=\frac{\hbar^{2}}{2m}+(C^{\tau}_{0}+C^{\tau}_{1})\rho\;\;\;\;\; \textrm{(PNM).} \end{equation} \tag{ 7 }$$

These coupling constants are C^τ₀ = 55.18 MeV fm⁵ and C^τ₁ = 1.23 MeV fm⁵ for the new SLy5* and C^τ₀ = 56.25 MeV fm⁵ and C^τ₁ = 23.95 MeV fm⁵ for the original SLy5. We observe that, while the two isoscalar coupling constants are very similar, the isovector one is very different. This, according to equation (6), explains the difference between the effective masses predicted by SLy5 and SLy5* as seen in figure 3. The fact that SLy5* gives a similar effective mass for the neutron in SNM and PNM is due to the very small value of C^τ₁ obtained with this interaction.

4.2. Binding energies

To further test the quality of the new SLy5* interaction, we calculated the binding energies of a set of semi-magic nuclei using the Hartree–Fock–Bogoliubov method. For the particle–hole channel we use the two functionals SLy5 and SLy5* while for the particle–particle channel we used a contact pairing interaction of the surface type [17].

With such a choice, we obtain an average pairing gap [18] Δ_aver = 1.201 MeV for SLy5 and Δ_aver = 1.164 MeV for SLy5* in ¹²⁰Sn. We compare the resulting binding energies with the experimental ones given by Audi et al [19]. We observe that there is no significant qualitative difference between the two functionals for this set of nuclei. The binding energies predicted with SLy5* are on average slightly smaller than those with SLy5 but both functionals compare with the experimental data with the same average deviation (figure 4).

Zoom In Zoom Out Reset image size

4.3. Charge radii

In figure 5, the charge radii obtained with the two functionals are compared with experimental results taken from the UNEDF webpage⁴. Once again, we observe that SLy5 and SLy5* show the same behavior, and the agreement between theory and experiment is of the same level.

Zoom In Zoom Out Reset image size

4.4. Stability against spin polarization

We have shown that, with the inclusion of the LR method into the fitting protocol, it is possible to build functionals that do not present instabilities. As a first example, we used this improved protocol to cure the spin instability present in SLy5. The resulting functional has a comparable level of accuracy as the original one for the prediction of relevant observables as binding energies and charge radii for the set of semi-magic nuclei considered here; moreover, most of the relevant infinite nuclear matter properties are unchanged; only the Thomas–Reiche–Kuhn enhancement factor is significantly larger.

Such a method clearly demonstrates that it is possible to use the LR method for additional constraints to build a functional with no finite-size instabilities in infinite nuclear matter and nuclei.

This work was supported by the Agence Nationale de la Recerche under grant number ANR BLANC 0407 'NESQ'. VH acknowledges financial support from the F.R.S.-FNRS (Belgium). The authors also thank T Duguet, P-H Heenen, M Bender for inspiring and interesting discussions.