Accurate relativistic density functional for exchange energy of atomic nuclei

The inclusion of nucleonic exchange energy has been a long-standing challenge for the relativistic density functional theory (RDFT) in nuclear physics. We propose an orbital-dependent relativistic Kohn-Sham density functional theory to incorporate the exchange energy with local Lorentz scalar and vector potentials. The relativistic optimized effective potential equations for the local exchange potentials are derived and solved efficiently. The obtained binding energies and charge radii for nuclei are benchmarked with the results given by the traditional relativistic Hartree-Fock approach, which involves complicated nonlocal potentials. It demonstrates that the present framework is not only accurate but also efficient.

Solving quantum mechanical many-body problems plays an essential role in many fields, such as materials science, condensed matter physics, nuclear physics, etc. Density functional theory (DFT) is one of the most successful tools, and no other method achieves comparable accuracy at the same computational costs. The key task for DFT is to build the a priori unknown energy density functional (EDF), whose existence is proved by the Hohenberg-Kohn theorem [1] but the exact form is always extremely difficult to determine and must be built with approximations. One of the most popular ways is provided by the Kohn-Sham DFT, where the EDF of an interacting system is constructed by introducing an auxiliary noninteracting system moving in a local potential V KS that gives the same ground-state density [2]. The EDF for the kinetic energy, the external potential energy, and the Hartree energy can be built straightforwardly, and the so-called exchangecorrelation energy functional for the remaining parts is to be addressed.
In nuclear physics, the DFT has been used with great success in describing various phenomena of nuclei throughout the nuclear chart [3,4]. The relativistic density functional theory (RDFT) [5,6,7,8] is of particular interest since it exploits basic properties of QCD at low energies, in particular, symmetries and the separation of scales [9]. By taking into account the Lorentz symmetry, the relativistic density functionals provide an efficient description of nuclei with a delicate interplay between the large Lorentz scalar and vector potentials of the order of hundred MeV, and explain the large spin-orbit splittings and the nuclear magnetic potential in a consistent way [10]. Over the past decades, a large variety of nuclear phenomena have been described successfully with the RDFT, ranging from infinite nuclear matter to spherical and deformed finite nuclei [5,6,7], from ground states to collective rotational and vibrational excitations [8,11], from static to dynamic properties [8,12,13,14].
In these studies, the underlying functionals contain only the kinetic and Hartree energies. The exchange energies are not taken into account, and the parameters in these functionals are phenomenological and cannot be directly related to realistic nucleon-nucleon interactions. Such a theoretical framework is relatively simple and the computational costs are low, but it brings following drawbacks: (1) the nuclear shell structures and their evolutions are not well reproduced due to the missing of the nucleonnucleon tensor interactions [15,16]; (2) an additional adjustment is needed to describe the charge-exchange spinflip excitations, such as Gamow-Teller and spin-dipole resonances [17,18]. Therefore, it is highly desirable to build a relativistic density functional including the exchange energies within the Kohn-Sham scheme. In fact, for Coulombic systems, it often allows the electronic EDFs to depend explicitly on the single-particle orbitals of the Kohn-Sham DFT [19,20], which are motivated by both practical and formal inadequacies in the conventional functionals, such as the presence of the self-interaction, the absence of a derivative discontinuity, etc. Such orbital-dependent functionals not only include exchange terms, but also provide a clear path toward an ab initio DFT. Indeed, the total energy can be generally written into orbital-dependent functionals in the many-body perturbation theory, which is directly connected to the microscopic Hamiltonian.
For nuclear systems, the relativistic many-body perturbation theory has achieved great successes in the past decade. In its simplest form, the relativistic Hartree-Fock (RHF) calculations reproduce successfully the shell structure evolutions [21,22,23,24,25] and spin-isospin resonances [17,18]. Moreover, in its resummed form, the relativistic Brueckner-Hartree-Fock (RBHF) theory has a direct connection to the realistic nucleon-nucleon interactions and can nicely reproduce the binding energies and charge radii of finite nuclei [26]. However, one has to introduce nonlocal potentials in RHF theories, so the theoretical framework is much more involved and the computational costs are extremely heavy [27,28]. With contact interactions, the RHF framework can be significantly simplified by representing the exchange terms as Hartree terms via the Fierz transformation [29,30], while the exchange of light π mesons cannot be considered since it is associated with the long-distance dynamics.
In this Letter, we construct, for the first time, an orbitaldependent relativistic density functional for the nuclear RHF energy, in which the single-particle RHF orbitals are replaced with the Kohn-Sham orbitals. In contrast to the complicated nonlocal potentials involved in the RHF approach, the present orbital-dependent RDFT works in the Kohn-Sham scheme and utilizes fully local relativistic Kohn-Sham (RKS) potentials. The RKS potentials are determined by the relativistic optimized effective potential (ROEP) method, which is implemented for nuclei for the first time in this work, while it has been widely used in Coulombic systems [31,32,33,34]. Unlike the Coulombic systems, nuclei are self-bound systems with both spin and isospin degrees of freedom, so the nuclear ROEP equations are strongly coupled in the Lorentz scalar and vector channels. Note that the nonrelativistic orbital-dependent DFT has been tested for artificial neutron drops based on the simple Minnesota nucleon-nucleon interaction with only central forces [35]. Nevertheless, the present density functional is fully relativistic and is based on a well-defined effective Lagrangian including nucleons, mesons, and photons. It can be used to describe real nuclei and we successfully benchmark our results against the conventional RHF results for nuclei from light to heavy.
The starting point of the present orbital-dependent RDFT is the RHF energy derived from an effective Lagrangian where the nucleons interact with σ, ω, and ρ mesons as well as the photon. Following the RHF theory [36,37], the RHF energy can be written as where M is the nucleon mass and the two-body interaction V (1, 2) includes the following meson-and photon-nucleon interactions: with D φ the propagator for meson and photon fields, g φ the coupling constants, e the charge unit, and τ the isospin Pauli matrices. The three terms in Eq. (1) correspond to the kinetic energy T , the Hartree energy E H , and the exchange energy E x , respectively. In contrast to the RHF theory, the occupied single-particle states |a and |b in the Kohn-Sham theory should be interpreted as the Kohn-Sham orbitals (defined in Eq. (5)). As a result, T and E x are orbitaldependent functionals, while E H can be written explicitly as a functional of the scalar density ρ s and the vector currents j µ that are defined with the Kohn-Sham orbitals in the coordinate space by, a∈τφ a (r)γ µ ϕ a (r).
Here, the subscript τ is used to distinguish neutron and proton.
According to the Kohn-Sham DFT, the RKS orbitals should be calculated by solving the RKS equation, which is essentially a Dirac equation, in which S τ and V µ τ are the RKS potentials defined via the variation of the energy functional, Each potential is local and is decomposed into the Hartree and exchange parts. Since the Hartree energy E H is an explicit functional of density and currents, one can readily obtain the Hartree potentials S H,τ (r) and V µ H,τ (r) by variation. However, it is nontrivial to obtain the exchange potentials S x,τ (r) and V µ x,τ (r) because the exchange energy E x explicitly depends on the Kohn-Sham orbitals rather than the density or currents.
The general procedure to get the exchange RKS potentials is to solve the ROEP equations, which are derived here via the chain rule of functional differentiation and the details can be seen in the Supplemental Material [38]. The nuclear ROEP equations consist of a pair of coupled equations for the scalar and vector potentials, occ.
Here,ξ a are the first-order changes in the RKS orbital ϕ a when the potential S x,τ + γ µ V µ x,τ in the RKS equation (5) is replaced with the orbital-specific potential W a (r) = 1 ϕa δEx δϕa . In this sense, the ROEP equations (8a) and (8b) are consequences of the vanishing first-order corrections for the density and currents. By introducing the first-order correction ζ a of the eigenvalue ε a , the ROEP equations can also be written equivalently as in which T es1 and T µ es2 are the energy-shift terms T os1 and T µ os2 are the orbital-shift terms Since ζ a andξ a correspond to the first-order perturbation of the energy ε a and the RKS orbital, respectively, the energy-shift and orbital-shift terms should be smaller than the W a terms on the righthand sides of Eqs. (9) and (10). As a demonstration, the present framework is applied to the spherical nuclei 16 O, 40 Ca, 48 Ca, 132 Sn, and 208 Pb. The calculated results are benchmarked with the RHF results with the widely used effective interaction PKO2 [21]. We start with an initial guess for the RKS potentials, and solve the RKS equation by expanding the RKS orbitals in a set of spherical Dirac Woods-Saxon (DWS) basis [39], which is constructed in a radial box with the size R box = 16 fm and the mesh ∆r = 0.1 fm. The RKS equations provides the RKS orbitals, which are used to calculate the scalar and vector densities, and in turn, the Hartree potentials. Then, the ROEP equations are solved to obtain the exchange parts of the RKS potentials. Note that the orbit-shift terms are always neglected here since they are quite small compared to other terms. This is consistent with the widely used relativistic Krieger-Li-Iafrate (RKLI) approximation [40,41] for Coulombic systems, while the validity of this approximation is confirmed for nuclear systems in the present work. Updating the Hartree and exchange potentials, the RKS equation is solved iteratively until convergence is achieved. In Fig. 1, our results for the total energy per nucleon as a function of the iteration steps for 16 O, 132 Sn, and 208 Pb are compared with the corresponding RHF results. During the iteration, the calculated energies with both methods smoothly converge to almost the same value with a similar convergence pattern. This clearly demonstrates the validity of the present orbital-dependent RDFT framework.
Although the iteration numbers to achieve convergence for the present RKS calculations are almost the same as the RHF calculations, the total computing time is much less; only about one third of the computing time of the RHF calculations, as can be seen in the inset clocks. Note that the RHF equation is a complex integro-differential equation with nonlocal exchange potentials, while all the potentials in the present RKS equation are local. As a result, it is not necessary to calculate the matrix elements of the nonlocal potential on the DWS basis in the solution of the RKS equation. This is the main reason of the computational merits, although this advantage comes at the price of the additional solution of the ROEP equations. In the present work, the ROEP equations are solved by neglecting the insignificant orbit-shift terms, but a preliminary full solution of the ROEP equations with the orbital-shift terms shows that the computing time does not change significantly. Note that a high accuracy full solution of the ROEP equations is quite a tricky problem, which has been discussed extensively in the field of Coulombic calculations [34,42]. Works along this direction are in progress. circles, and open triangles denote the RHF results, the RKS results, and the RKS results neglecting the energy-shift terms in the ROEP equations, respectively. Differences in the total energies (b) and the charge radii (c) between the RKS and RHF results with and without the energy-shift term are also shown. Note that the microscopic center-of-mass correction energies [43,44,45] are included here.
The influence of the orbit-shift and energy-shift terms on the calculated total energies and charge radii is shown in detail in Fig. 2. The ground-state binding energies and charge radii of 16 O, 40 Ca, 48 Ca, 132 Sn, and 208 Pb are calculated within the present RKS framework, and are compared with the RHF results. Note that the microscopic center-of-mass correction [43,44,45] is considered in both RKS and RHF calculations after the iteration converges. As can be seen in Fig. 2(a), the RHF results are well reproduced by the present RKS calculations where the orbit-shift terms are neglected. By further neglecting the energy-shift terms, the RKS results are still in a reasonable agreement with the RHF results, and visible deviations can be seen only for heavy nuclei such as 132 Sn and 208 Pb. This demonstrates that the W a terms in the ROEP equations alone provide a good estimation for the local exchange potentials.
A more detailed comparison is made in Figs. 2(b) and 2(c) by depicting the differences in the total energies and charge radii between the RKS and RHF results with and without the energy-shift term. Although both the RKS and RHF calculations aim to minimize the RHF energy, the former should be regarded as a constrained optimization of the total energy by requiring local auxiliary potentials only. Therefore, the ground-state energies obtained via RKS are expected to be higher than those of RHF, and this is consistent with our calculations as shown in Fig. 2(b). The small differences between the RKS and RHF results demonstrate that the orbit-shift terms can be safely neglected for nuclear systems. The energy-shift terms, however, are relatively more important in particular for heavy nuclei. For example, the energy difference between the RKS and RHF results can reach to 5.7 MeV for 208 Pb. The main difference between the present RKS and the RHF approaches is that the latter involves nonlocal potentials, whereas the former employs only local Kohn-Sham potentials. In Fig. 3, the nucleon and anti-nucleon RKS potentials for 16 O are depicted in comparison with the RHF mean potentials. The Hartree potentials are local for both RKS and RHF approaches, and they nicely agree with each other. The magnitudes for the Hartree potentials are around several tens MeV in the Fermi sea and several hundreds MeV in the Dirac sea, which is associated with the large spin-orbit splittings in nuclei. The exchange potentials in the RHF method are nonlocal and cannot be directly compared to local potentials. In the present RKS framework, however, the exchange potentials are fully local. It is seen that the magnitudes of the exchange potentials are smaller than those of the Hartree potentials in both the Fermi sea and Dirac sea. In particular, in the Fermi sea, the obtained exchange potential is slightly repulsive in contrast to the attractive Hartree potential. As a result, they provide opposite contributions to the total energy. It is of importance to have the correct long-range behavior of the Coulomb potential for protons, i.e., ∼ (Z − 1)e 2 /r (Z is the proton number). This is in fact not a trivial condition and the Hartree potential decays as ∼ Ze 2 /r due to the self-interaction. Figure 4 shows the proton potentials of 16 O at large radial distances given by the RKS calculation. By including the exchange energy that removes the self-interaction, we are pleased to find that the RKS calculation nicely reproduces the correct long-range behavior.
In summary, an orbital-dependent relativistic energy density functional including the nucleonic exchange energy has been constructed for the first time within the relativistic Kohn-Sham framework for the nuclear RHF energy. The single-particle RHF orbitals are replaced with the Kohn-Sham orbitals, which are obtained by solving the relativistic Kohn-Sham equation with local Lorentz scalar and vector potentials. Although it requires an additional solution of the ROEP equations to obtain the local exchange potentials, the present framework is superior to the traditional RHF approach which involves complex nonlocal potentials. Benchmark calculations have been performed for the ground-state energies and charge radii of nuclei from light to heavy, and it demonstrates that the proposed framework is not only highly accurate but also efficient.
Since the exchange potentials are local, this new framework could be straightforwardly extended to nuclei with arbitrary deformation by solving the RKS equations on three-dimensional lattice [46,47]. Moreover, it provides a promising way to construct an ab initio relativistic energy density functional for atomic nuclei based on realistic nucleon-nucleon interactions, such as the Bonn potentials [48] and the relativistic chiral interactions [49,50], softened with the Brueckner G-matrix approach [51,52] and/or modern renormalization group methods [53].