Ab initio construction of the energy density functional for electron systems with the functional-renormalization-group-aided density functional theory

We show an $\textit{ab initio}$ construction of the energy density functional (EDF) for electron systems using the functional renormalization group. The correlation energies of the homogeneous electron gas given in our framework reproduce the exact behavior at high density and agree with the Monte-Carlo data in a wide range of densities. Our analytic technique enables us to get the correlation energies efficiently for various densities, which realizes the determination of EDF in the local density approximation (LDA) without any fitting for physically relevant densities. Applied to the Kohn-Sham calculation for the noble gas atoms, our EDF shows comparable results to those of other conventional ones in LDA.

as a practical method. In particular, as systems for which DFT is frequently employed, the electron systems are one of the most important targets.
The aim of this Letter is the microscopic derivation of the EDF E [ρ] for the spin-unpolarized electron systems with the aid of the FRG-DFT. To our best knowledge, this work is the first numerical application of the FRG-DFT to (3 + 1)-dimensional systems. As a first step of the microscopic construction of the EDF, we consider the exchange-correlation part in the local density approximation (LDA) and aim at the construction of the correlation part. To this end, we apply the FRG-DFT to the (3 + 1)-dimensional homogeneous electron gas (3DHEG), derive the expression for the correlation energy per particle ε corr by solving the flow equation analytically with employing the second-order vertex expansion, and obtain ε corr (r s ) as a function of the Wigner-Seitz radius r s = [3/ (4πρ)] 1/3 . Our ε corr (r s ) reproduces the exact behavior at the high-density limit given by the result of Gell-Mann-Brueckner resummation ε GB corr (r s ) [27] and agrees with the results of the diffusion Monte-Carlo (DMC) calculations [28][29][30] in a wide range of densities.
A concern about the application to (3 + 1)-dimensional systems may be that the coordinate or momentum integrals in the FRG-DFT calculation become time consuming. For 3DHEG, however, we find that the dimension of the integrals can be drastically reduced with analytic techniques, which enables us to obtain ε corr (r s ) densely enough to determine the EDF without fitting for physically relevant densities. This is in contrast to other conventional LDA EDFs, most of which are determined based on empirical choice of the fitting function for a few DMC data points.
Furthermore, applying the KS calculation of the ground states of the noble gas atoms, we demonstrate that the EDF constructed from our FRG-DFT data shows comparable results to other conventional LDA EDFs.
In this Letter, Hartree atomic units are employed. FRG-DFT . We briefly summarize the formalism of FRG-DFT and our analytic results in the 3DHEG, where electrons and background ions neutralizing the system interact to each other via the two-body Coulomb interaction U (x) = 1/ |x|. Following Refs. [6,7], we consider the evolution when the inter-particle interaction is gradually turned on. Let us employ a parametrized two-body interaction U λ (x) with the evolution parameter λ running from λ = 0 to 1 and set U λ=0 (x) = 0 and U λ=1 (x) = U (x) to describe the evolution from the free to the fully interacting systems. Not only the electronelectron but also the electron-ion and ion-ion interactions are substituted by U λ (x) so as to keep the system neutral and avoid divergence caused by the Hartree energy during the evolution [24].
The key quantity of the FRG-DFT is the effective action Γ λ [ρ] for density ρ. To define it, we start from the action depending on λ in the imaginary-time formalism: Here, we have introduced X = (τ, x) and X = dτ dx with imaginary time τ and spatial coordi- , ψ ↓ (X)) standing for the electron field with spin ↑ and ↓, andρ ∆ (X) =ρ (X) − n e with the density fieldρ (X) = ψ † (X ε ) ψ (X) and n e = 3/ 4πr 3 s being the densities of electrons and back-ground ions. The second term includes the electron-electron, electron-ion, and ion-ion interaction terms. We have also introduced X ε = (τ + ε, x) with an infinitesimal ε > 0 so that the Hamiltonian corresponding to S λ becomes normal ordered [23]. Then, S λ defines Γ λ [ρ] as is the generating functional for density correlation functions and ρ(X) is an arbitrary density. A notable feature of Γ λ [ρ] is that it satisfies the variational principle and gives the ground-state energy and density [17,19] [7,20,21,23,24]: where ρ ∆ (X) = ρ (X) − n e and Γ Also, X ε is defined in the same manner as X ε but ε → 0 limit is taken after ε → 0 so that Γ (2)−1 λ [ρ] (X ε , X ) can be treated as the density correlation function [23]. The crucial point of Eq. (1) is that it is written in a closed form of Γ λ [ρ], which provides systematic schemes for the derivation of Γ λ [ρ].
Practically, the functional differential equation (1) needs to be converted to some numerically solvable equations. Here, we introduce the vertex expansion [6,7]: The functional Taylor expansion around a homogeneous density ρ (X) = n e is applied to Eq. (1), which yields a hierarchy of differential equations for density correlation functions [21,[23][24][25]. We consider the expansion up to the second order and truncate higher-order terms. The equations up to the second order in the momentum-space representation read ∂ λG whereŨ λ (p) is the Fourier transform of U λ (x) and ε gs,λ = lim β→∞ Γ λ [n e ] / (βN ) with N = n e dx being the total particle number is the ground-state energy per particle. Here, we have introduced P = (ω, p) and 3 with the Matsubara frequency ω and the spatial momentum p, andG (n) λ (P 1 , . . . , P n−1 ) being the connected density correlation function. Since C λ (P ) is composed of higherorder correlation functions, an approximation for it is required. Here, we employ the approximation C λ (P ) ≈ C λ=0 (P ), with which Eqs. (2) and (3) can be solve analytically when U λ (x) is chosen as U λ (x) = λU λ=1 (x) [24].
Extracting ε corr from ε gs,λ=1 , we obtain which plays the central role in our construction of the EDF. Here, we have introduced f (x, y) = cosh y + (x/y) sinh y, A P :=Ũ (p)G (2) λ=0 (P ), and B P := Ũ (p) C λ=0 (P ) 1/2 , in whichŨ (p) = 4π/p 2 is the Coulomb interaction in the momentum representation and C λ=0 (P ) is evaluated from the connected density correlation function in the free system: with the symmetric group S n−1 of order n − 1, the two-point propagator of free fermionsG , and the Fermi momentum p F = (9π/4) 1/3 /r s . Before ending the summary of the formalism, we comment on the behavior of Eq. (5) at the high-density limit r s → 0: Through the use of the scaling behavior with the dimensionless momentum P = r 2 s ω, r s p , the expansion of Eq. (5) with respect to r s is obtained. From the expansion, one finds that the exact behavior ε GB corr (r s ) is reproduced at r s → 0.
Reduction of dimension of multi-integral . A difficulty in the three-dimensional system may be that the numerical evaluation of multi-integrals with respect to momenta is too costly to get ε corr (r s ) for various r s . We find that this can be circumvented since the dimension of the integral in Eq. (4) can be drastically reduced in an analytic manner.
Using the expression forG (n) λ=0 (P 1 , . . . , P n−1 ) given in the sentence below Eq. (5) and performing the frequency integral, Eq. (4) becomes where and θ p = θ (−ξ (p)) with the Heaviside step function θ (x). By usingŨ (p) = dx e ip·x / |x| and employing the cylindrical coordinates with choosing the direction of p as the direction of the longitudinal axis (z axis), Eq. (6) is rewritten as follows: Here, we have introduced P r (p z ) = p 2 F − p 2 z 1/2 , D (ω, p, p z ) = iω + pp z + p 2 /2 −1 , and where J 1 (x) and K 0 (x) are the Bessel function of the first kind and the modified Bessel function of the second kind, respectively. Then, the integral in I (a, b, c) can be performed analytically [31]: where Finally, Eq. (7) together with Eq. (8) shows that only a double integral is required for the calculation of C λ=0 (P ). Needless to say, the isotropy reduces the dimension of the integral in Eq. (5).
Results of the correlation energy. The reduction of the dimension of the integral saves the time for the numerical calculation and enables one to obtain ε corr for many r s . The calculation was carried out on 65536 grid points with the logarithmic mesh in r s ∈ 10 −6 a.u., 100 a.u. . Figure 1 shows ε corr (r s ) obtained by the FRG-DFT together with ε GB corr (r s ) and the DMC results [28][29][30]. As expected from our analytical discussion, one can see the FRG-DFT result reproduces the exact behavior at the high-density limit given by ε GB corr (r s ). The FRG-DFT result is also consistent with the DMC results in a wide range of r s , and in particular the agreement becomes better as the density increases: The deviations from the results in Ref. [28] are about 2.0 % at r s = 1 a.u., 6.7 % at r s = 50 a.u., and 17 % at r s = 100 a.u. Now we have shown that ε corr (r s ) can be obtained in the framework of the FRG-DFT and becomes accurate as the density increases, which is one of the main results of this Letter. In the remaining part, we attempt a construction of the EDF by use of our ε corr (r s ).
Construction of the energy density functional . Since ε corr are obtained very densely for various r s in our The DMC results are obtained by subtracting the kinetic and exchange energies [32,33] from the total energies given in Refs. [28][29][30].
scheme, we can construct the LDA EDF E LDA corr [ρ] = dx ρ (x) ε corr (r s (ρ (x))) without any fitting for physically relevant densities. This is in sharp contrast to other conventional EDFs such as VWN [34], PZ81 [35], and PW92 [36], which are determined by fitting the few DMC data obtained by Ceperley and Alder [28].
Our functional, which is referred to as the FRGnumerical-table functional (FRG-NT), is constructed as follows: In r s ∈ 10 −6 a.u., 100 a.u. , ε corr (r s ) are determined by the interpolation of the FRG-DFT data. For simplicity, we employ the linear interpolation; the results hardly depend on the choice of the interpolation function. In r s < 10 −6 a.u., ε corr (r s ) is substituted by ε GB corr (r s ). The FRG-DFT data are extrapolated to r s ≥ 100 a.u. by a fitting function ε corr (r s ) = γ/ 1 + β 1 √ r s + β 2 r s [35].
The fitting parameters are chosen to be γ = 0.0378052, β 1 = −0.801035, and β 2 = −0.0306778, which are obtained by fitting the data in 95 a.u. < r s < 100 a.u. A remark is in order here: The evaluation of ε corr (r s ) = dε corr (r s ) /dr s appearing in the KS potential δE LDA corr [ρ] /δρ (x) = [ε corr (r s ) − (r s /3) ε corr (r s )] rs=rs(ρ(x)) with the numerical differentiation may cause numerical errors. We evade this by performing the analytic differentiation of Eq. (5): where g (x, y) = x+(x cosh y + y sinh y) /f (x, y). To perform the differentiation, we have usedG  dent of r s . We calculate ε corr (r s ) on the same grid for r s as ε corr (r s ) and determine ε corr (r s ) for arbitrary r s in the same manner as ε corr (r s ). Additionally, we prepare a functional, which we name FRG-PZ, by fitting the FRG-DFT data with the same function as PZ81 ε corr (r s ) = A ln r s + B + Cr s ln r s + Dr s r s < 1 a.u., γ/ 1 + β 1 √ r s + β 2 r s r s ≥ 1 a.u., (10) for the purpose of comparing our functionals to PZ81 and discussing the origin of the deviation between EDFs. Here, A = 0.0311 and B = −0.0480 reproduce ε GB corr (r s ) at r s → 0 . The remaining parameters C, D, γ, β 1 , and β 2 are related to each other through the continuum conditions for ε corr (r s ) and ε corr (r s ) at r s = 1 a.u: Table I lists the values of the parameters obtained by the fitting with the conditions Eqs. (11a) and (11b). Benchmark test of the functionals. We apply our EDF to the KS calculation of the ground-state energies of noble gas atoms and compare with other conventional EDFs such as VWN [34], PZ81 [35], PW92 [36], Chachiyo [37], revChachiyo [38], and GGA-PBE [39]. The numerical calculation was carried out by use of ADPACK [40]. Figure 2 shows the ground-state energies E gs of Ne, Ar, Kr, Xe, and Rn atoms obtained by each EDF as ratios to the results of PZ81 E gs PZ81 . One can see that the functionals constructed from the FRG-DFT show comparable results to those of other LDA EDFs for every atoms. On the other hand, the GGA-PBE results are quite different from those of the LDA EDFs. This suggests that the results of our functional reside near that of the exact LDA EDFs, i.e., the LDA EDF constructed from the exact correlation energy per particle ε exact corr (r s ), as much as other conventional LDA EDFs, and the differences between our functional and other LDA EDFs are insignificant for the accuracy in comparison with the effect of the ignorance of the gradient.
Origins of the difference among functionals. To further understand the difference between LDA EDFs, we inves- tigate the difference of the energy with fixing the density inspired by the notion of the functional-driven error [41]. Figure 3 shows ∆E gs gs (x) being the ground-state density obtained by PZ81. One can see that the deviations among EDFs are comparable even at the same density. In the case of LDA, this deviation originates from ε corr (r s ), the error of which is expected to be attributed to two parts: the reference-driven error and the fitting-driven error, i.e. the errors caused by the choice of the reference data and the fitting functions, respectively. By recasting the difference between EDFs in terms of these two errors, further understanding of the origin of ∆E gs F shown in Fig. 3 will be obtained.
We roughly assume that ε i corr (r s ) standing for ε corr (r s ) used for functional i (= VWN, PZ81, PW92, FRG-NT, FRG-PZ) is written as  Figure 4 shows ∆ε i corr,fit (r s ) and ∆ε FRG corr,ref (r s ) − ∆ε DMC corr,ref (r s ). In r s 10 a.u., the use of fitting affects the value of ε corr (r s ) more than the choice of the reference data, while these quantities have comparable magnitude in r s 10 a.u. This may explain why the comparable results for each functional are obtained in Fig. 3 since r s 10 a.u. is physically relevant for atoms.
Conclusion. We presented an ab initio construction of the energy density functional (EDF) for threedimensional electron systems using the functionalrenormalization-group-aided density functional theory (FRG-DFT). The derived correlation energies of the homogeneous electron gas agree with the Monte-Carlo results in a wide range of densities reproducing the exact behavior given by the Gell-Mann-Brueckner resummation at the high-density limit. Using the FRG-DFT data obtained densely for various densities, we construct the EDF in the local density approximation (LDA) without using any fitting function for physically relevant densities. Applied to the KS calculation of the ground-state energies of the noble gas atoms, our functional shows comparable results to other conventional ones in LDA. Our results show that FRG-DFT can become a practical method contributing to the non-empirical construction of EDFs of realistic quantum many-body systems.
Although we have focused on the case of LDA in this Letter, our formalism is also applicable for the construction of EDFs incorporating the effect of the gradient of density. There are some technical ideas to realize the inclusion of gradient effects, such as the use of the weighted density approximation or the derivative expansion, which has been developed in the context of FRG. Based on these ideas, we believe that the construction of EDFs beyond LDA without any empirical parameter is achievable. Our formalism and procedure can also be naturally extended to the case of constructing EDFs in the local spin density approximation, which will be presented in a forthcoming paper.