Kinetic and Exchange Energy Densities near the Nucleus

We investigate the behavior of the kinetic and the exchange energy densities near the nuclear cusp of atomic systems. Considering hydrogenic orbitals, we derive analytical expressions near the nucleus, for single shells, as well as in the semiclassical limit of large non-relativistic neutral atoms. We show that a model based on the helium iso-electronic series is very accurate, as also confirmed by numerical calculations on real atoms up to two thousands electrons. Based on this model, we propose non-local density-dependent ingredients that are suitable for the description of the kinetic and exchange energy densities in the region close to the nucleus. These non-local ingredients are invariant under the uniform scaling of the density, and they can be used in the construction of non-local exchange-correlation and kinetic functionals.


I. INTRODUCTION
Kohn-Sham (KS) density functional theory (DFT) 1-3 can be considered the most used method in electronic calculations of quantum chemistry and condensed matter physics. Its practical implementation is based on approximations of the exchange-correlation (XC) energy (E xc ), which is a subject of intense research 4-6 . Moreover, subsystem DFT 7-9 and orbital-free DFT 10-13 need the use of kinetic energy (KE) functional approximations.
Concerning XC functionals, the simplest ones beyond the local density approximation (LDA) are those based on the generalized gradient approximation (GGA) [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] , which are constructed using the electron density (ρ) and its gradient (∇ρ). Meta-generalized gradient approximations (meta-GGAs) [31][32][33][34][35][36][37][38][39][40][41][42] are the most sophisticated semilocal functionals, incorporating important exact conditions and having an improved overall accuracy with respect to the GGA functionals. The meta-GGA functionals use as an additional ingredient to the GGA ones the Kohn-Sham positive KE density: (the total KE being T = τ KS d 3 r). In Equation (1), f = 2 for closed-shell systems (as the ones considered in this work), and the summation is over all occupied orbitals; atomic units, i.e., e 2 = = m e = 1, are used throughout. The quantity τ KS enters in the expansion of the angle-averaged exact exchange hole 4, 43 , being thus a natural and important tool in the construction of XC approximations. An important requirement for an accurate XC functional is a proper model for the exchange energy (XE) density 15,44 . The XE density is usually defined in terms of the exchange enhancement factor: where e HEG x = −(3/(4π))(3π 2 ρ) 1/3 ρ is the XE density of the homogeneous electron gas (HEG). We recall that the XE density is not uniquely defined (being up to a gauge transformation), but its underlining hole must be realistic and close to the exact one, which is an observable [45][46][47][48] . In this work, we will use as reference the definition of the conventional exact-exchange density. At the GGA level, the total exchange energy is usually expressed as: i.e., with the exchange enhancement factor being a function of the reduced gradient s = |∇ρ|/[2k F ρ]. Here, k F = (3π 2 ρ) 1/3 is the local Fermi wavevector. At the meta-GGA level, we have: i.e., F x is also a function of reduced Laplacian q = ∇ 2 ρ/[4k 2 F ρ] and/or of τ KS . At the nucleus of the helium isoelectronic series, it has been shown 49 that F x = 1 3 ( 4π 2 3 ) 1/3 ≈ 0.787. However, in popular semilocal exchange functionals (GGAs and meta-GGAs), F x ≥ 1: thus, these functionals cannot be realistic near the nucleus, where there is an important de-enhancement. The nuclear region can be identified using the usual semilocal ingredients. For example, the reduced gradient s behaves at the nucleus of the helium isoelectronic series as s = 1/(6π) 1/3 ≈ 0.376, while the reduced Laplacian q diverges to −∞. This issue has been considered by Tao 49 , who constructed an exchange functional with the correct XE density at the nucleus, using the single inhomogeneity parameter proposed by Becke 43 Q B = 1 − τ τ HEG + 5s 2 3 + 10 3 q with τ HEG = (3/10)(3π 2 ) 2/3 ρ 5/3 . Further improvements on the development of Laplacian-dependent exchange functionals have been found by Cancio et al. 50 .
Similar shortcomings as for the XE density near the nucleus affect also many KE functionals at the GGA level [51][52][53][54][55] or at the Laplacian level 56 ; for a recent review of semilocal functionals, see 57 . The KE density is usually defined in terms of the KE enhancement factor: so that the total kinetic energy is: The von Weizsäcker (VW) kinetic energy density: (i.e., F W s = (5/3)s 2 ) is expected to be very accurate at the nucleus 2,58-60 . At the nucleus of the helium isoelectronic series, F s ≈ F W s ≈ 0.2353, i.e., there is an even more pronounced de-enhancement than in the exchange case. On the other hand, most of the GGA KE functionals have F s ≥ 1, while few semilocal KE functionals recover the VW at the nucleus 59,61,62 .
However, we recently pointed out that the VW functional does not have the correct behavior at the nucleus 63 . It was proven that the KE density at the nuclear cusp behaves as 63 : where ρ s and ρ p are the densities of s-type and p-type shells, respectively. Thus, also p-shells contribute to the KE density at the nuclear cusp 63,64 . The second term on the right-hand-side of Equation (8) has been evaluated for real atoms, and its contribution in the semiclassical limit of a neutral atom with an infinite number of electrons reaches 12% of the total KE density 63 . Understanding the physical phenomena at the nucleus can thus boost the development of more accurate XC and KE approximations. We recall that the nucleus region contains an important part of the total kinetic and exchange energies, and thus, small modifications of the KE and XC enhancement factors can bring significant variations to total energies.
In this paper, we will consider different aspects of density functionals in the nuclear regions: (i) we will describe the difference between the exact and VW kinetic energy densities, in a region near the nucleus, extending the derivation of 63 , where only the nuclear cusp was considered; (ii) we will present an approach based on the helium iso-electronic series, which correctly describes the nuclear region, for small atoms up to the semiclassical limit of large non-relativistic neutral atoms; (iii) we will propose novel non-local density ingredients for the conventional exchange and kinetic energy densities at the nuclear region.
The paper is organized as follows: In Section II, we present a detailed analysis of the kinetic and exchange energies at the nucleus of spherical systems.
We investigate the hydrogenic shells, the ten-electron hydrogenic model and the asymptotic neutral atom with an infinite number of electrons, presenting the small-r expansions of various quantities of interest. We also demonstrate that the 1s-shell model (1SM) approach is remarkably accurate for both kinetic and exchange energies, in the case of real atoms. The simple 1SM cannot be described by any semilocal ingredient, having a significant amount of non-locality. Consequently, in Section III, we propose non-local ingredients for exchange and kinetic energies, near the nucleus. These approximations are invariant under the uniform scaling of the density, and they can be used in the construction of non-local functionals. Finally, in Section IV, we summarize the results.

II. KINETIC AND EXCHANGE ENERGY DENSITIES AT THE NUCLEAR CUSP IN SPHERICAL SYSTEMS
For a system in a central potential, the KS orbitals can be written as φ nlm (r) = R nl (r)Y lm (θ, φ), where R nl (r) are the normalized radial functions, Y lm (θ, φ) are spherical harmonics, n is the principal quantum number, l is the angular momentum and m is the azimuthal quantum number. The density of the shell nl is: where in the last equality, we used Unsöld's theorem for spherical harmonics. The Kohn-Sham KE density of the shell nl satisfies the relation 63,65,66 : which is valid at any radial distance r from the nucleus. Note that the above relations are valid for systems in which a single shell is occupied. For real systems, with many occupied shells, the total density and the total KE density are obtained summing over all shells, i.e.: Note that, instead, in general, τ W = nl τ W nl due to the non-linearity of the VW functional 63 .
We will also consider the exact-XE density (per volume), which is given by: where λ k l,l = (2l+1)(2l +1) 2k+1 | l0l 0|k0 | 2 , lml m |l m are Clebsch-Gordan coefficients, r < = min(r, r ) and r > = max(r, r ). In the relevant case of a system with only s-shells occupied, we have: The above formulas are completely general and apply to any electronic system with a central external potential, e.g., real atoms. Nevertheless, a very interesting special case is that of hydrogenic orbitals. In fact, in this case, all calculations are analytical, and explicit formulas can be obtained by any symbolic computer algebra system software. Moreover, we recall that while in real atoms the electrons far from the nucleus experience a screened nuclear charge so that the corresponding orbitals differ from the hydrogenic ones, for large atoms or very positive ions, this screening effect becomes vanishingly small, and the simple model of hydrogenic orbitals becomes exact 68 . This model system has been largely used in DFT [69][70][71][72] , is very important for semiclassical physics 69,73,74 and has been used as a main reference system in the APBE 20 and APBEk 55 GGA functionals.

A. Hydrogenic Shells
Contributions near the nucleus are given by s-type shells for the KE and the XE density and by p-type shells for the KE density only. Higher angular momenta do not contribute near the nucleus.
(i) For a filled s-shell (l = 0, for any principal quantum number n) with f = 2 electrons, we find (after some algebra) near the nucleus (of charge Z): and thus, Kato's theorem 75 : is satisfied for any n. Equations (16) and (17) show that at the nucleus, all s-electrons are important, even if the main contribution is given by the n = 1 term (i.e., the 1s-shell), due to the n 3 term at the denominator.
In case of exchange, we find that the XE density of the shell n0 is: Equation (20) shows that at the nucleus, only s electrons with very small n contribute, due to the n 5 term at the denominator. Expression (20) generalizes the one in 49 .
For the special case of n = 1 (i.e., the helium isoelectronic series with fixed hydrogenic orbitals), we obtain the following expressions for the kinetic and exchange enhancement factors: F HY D1s which can be seen as simple semilocal conditions at the nucleus; see also 49 . Equation (22) shows a significant de-enhancement (F x < 1) at the nucleus, which is not reproduced by conventional DFT functionals (all GGAs and most meta-GGAs).
(ii) For a p-shell (l = 1, for any principal quantum number n), we find near the nucleus (of charge Z): For a system with only the p-shell occupied, there is no cusp of the density, and Kato's formula cannot be applied. Interestingly and importantly, even if ρ n1 (0) = 0 at the nucleus, the kinetic energies are of the same order of magnitude (∼ Z 5 ) as in the s-shell case. From Equations (24) and (25), we find: for any principal quantum number, in agreement with Equation (8).
Similarly, the XE density at the nucleus, for the pshells, is: and thus, for the exchange case, the p-orbitals do not contribute at the nucleus.

B. Ten-Electron Hydrogenic Model
Now, let us consider a 10-electron hydrogenic atom, with electronic structure 1s 2 2s 2 2p 6 and nuclear charge Z (i.e., the Ne isoelectronic series with fixed hydrogenic orbitals). The calculations are analytical, and we obtain (after some algebra): Note that Equation (29) confirms the validity of Equation (8), and Equations (31) and (32) confirm that the VW approximation underestimates the exact result at the nuclear cusp and in a region close to the nucleus (the coefficients of Zr are different).
In the case of exchange, we obtain: Note that the two leading terms in Equations (33) and (34) are the same in the case of a Be isoelectronic series, since the p-shell contributes to the exchange only with power r 2 or higher.

C. The Asymptotic Neutral Atom with an Infinite Number of Electrons
Using Equations (11) and (12) and the Riemann ζfunction (i.e., ζ(s) = ∞ n=1 1/n s ), we obtain the following analytical expressions near the nucleus: As expected, the difference between the exact KE density and the VW one is larger than in the case of ten electrons: interestingly, the coefficients of Zr differ more significantly than the values in r = 0. Note that at the nuclear cusp of an atom with Z → ∞ electrons, the reduced gradient is s = 0.3534, i.e., smaller than for the helium isoelectronic series (s = 0.375).
For the exchange case, no analytic results could be found. In fact, even if only the s-type shell contributes to the nuclear cusp, Equation (15) involves double sums, and we could not find a closed form expression for the terms with n = n . We performed the double sums numerically up to n = n = 50 obtaining:

D. 1s-Shell Model
Comparing the results of the previous three sections, we note that the exact KE enhancement factors F s are very similar in all three cases (see Equations (21), (31) and (39)), with a maximum deviation of only 0.0014. This is not the case for F W s (deviation 20-times larger). This means that the KE enhancement factor near the nuclear cusp of an atom with an infinite number of electrons is almost equivalent to the one obtained from the 1s shell only.
We can thus estimate the KE density near the nuclear cusp via a simple procedure, which we call the 1s-shell model (1SM): for a given atom, we consider only the density form the 1s-shell and use it to compute the KE density (which then equals the VW one), i.e., For the neon hydrogenic isoelectronic series, the 1SM gives very accurate results: this, however, traces back to a subtle error cancellation between the 2s and 2p contributions.
It is important to underline that, despite the simple expression in Equation (42), the ρ 1s can hardly be described by any semilocal ingredient of the total density ρ: in other words, the mapping ρ → ρ 1s is highly nonlocal. We can conclude that F s ≈ 0.235 at the nuclear cusp will yield accurate results for all atoms, from He to the semiclassical limit, whereas F s ≈ F W s is not accurate. Concerning the exchange, comparing Equations (22), (34) and (41), we see that all of the exact exchange enhancement factors F x are similar to each other.
Thus, also for the exchange case, we can define the 1SM approach as: A similar approach has been used to compute the exchange potential at the nucleus 76 .
In the next section, we will show that the 1SM model will yield good accuracy for both the KE and XE density, also in the case of real atoms.

E. Real Atoms
We considered non-relativistic neutral noble atoms using self-consistent numerical Kohn-Sham exact exchange orbitals and densities 69,77 . The wavefunctions R nl are discretized on a semi-logarithmic numerical grid, using the Engel code 78,79 .
In Figure 1, we show the error on the kinetic enhancement factor F s − F exact s versus the scaled radial distance r/R near the nucleus of the Kr atom: we consider the VW functional, the Kato approximation (i.e., τ Kato = Z 2 ρ/2) and the 1SM method. The scaled radial distance is defined as the average distance of the 1s shell: where ρ 1s is the density of the 1s shell (and drρ 1s = f = 2). For real atoms, we find that R ≈ 3/(2Z). Figure 1 shows that: (i) the VW approximation is not exact at the nuclear cusp; (ii) the Kato expression is a very good model for the VW behavior, not only at the nuclear cusp, but also for r/R < 0.2 80 ; (iii) the 1SM is accurate at the nucleus. Figure 2 reports F s at the nuclear cusp, for neutral noble atoms with filled shells with the number of electrons in the range 2 ≤ Z ≤ 2022 63 , considering the exact F s , the VW functional and the 1SM approach. Note that the value F s (r = 0) has been extrapolated from the available numerical grid point closest to the nucleus. For the He atom, all approaches coincide (F s ≈ 0.3): the results differ from the F HY D1s s of Equation (21) due to the screening effects, which are largest for the He atom, but rapidly decrease with increasing Z. For larger atoms, the exact F s converges to 0.2367 (see Equation (39)), while F W s is much lower. On the other hand, the 1SM method nicely reproduces the exact results for almost all systems, converging to the HYD1s value (see Equation 21) for Z → ∞.
We now turn to the exchange case. In Figure 3, we show the exact exchange enhancement factor at the nu- cleus, the 1SM approach and the HY D1s value of Equation (22), for noble atoms (2 ≤ Z ≤ 2022). The 1SM approach yields very accurate results only for the smallest atoms (up to Ar). For the largest atoms, the differences are significantly larger than in the kinetic case: in fact, the XE density is more non-local than the KE one. The simple expression F HY D1s x (0) = 0.787 is also accurate and can be used for the construction of more realistic semilocal exchange functionals 49 .

III. NON-LOCAL APPROXIMATIONS FOR EXCHANGE AND KINETIC ENERGIES, AT AND NEAR THE NUCLEUS
As shown above, the kinetic and exchange energy densities are fully non-local near the nucleus, and thus, their behaviors cannot be well captured by semilocal ingredients.
In order to build non-local ingredients to describe better the features of the exchange and kinetic functionals, first we need to consider appropriate lengths. In Figure 4, we compare the VW and the Fermi wavelengths, defined by: The VW length is similar for all atoms that contain porbitals (Ne-Rn) and is slightly different for the He atom. Thus, λ W can distinguish between atoms that contain only s-orbitals (e.g., He) and the other atoms (e.g., Ne-Rn). For this reason, we chose this length for the kinetic case.
On the other hand, the Fermi wavelength near the nucleus is similar for all of the atoms (He-Rn) reported in the figure. We will use this length, combined with the other meta-GGA ingredient α = (τ KS − τ W )/τ HEG = F s − F W s , for the exchange case. In Figure 5, we report α which shows that α is not vanishing at the nucleus. Figure 5 shows that α has a monotonic behavior for the noble atoms series, starting from α = 0 in the case of the He atom, while for the Rn atom, it becomes quite close to the asymptotic limit. We will use this finding in order to construct a proper ingredient for exchange.

A. Kinetic Energy
In this section, we build a new non-local density ingredient suitable for the description of the KE density near the nucleus, i.e., not only at the nuclear cusp, but also in a region around the nucleus.
We start from the observation that the 1SM model of Equation (42) is almost exact; see Figure 2. Then, instead of considering the VW functional of the density of the 1s-shell, we consider a screened VW functional of the total density.
Thus, we propose the following β-family of ingredients: where a, b and β are positive constants. The non-local ingredient y W is invariant under the uniform scaling of the density (as s and q), and thus, it can be used in the construction of KE enhancement factor approximations. For simplicity, in this paper, we chose β = 1. Note that the exponential damping factor of Equation (48) uses the VW length of Equation (45).
The parameters a and b have been optimized for the Ne atom, in the region r/R ≤ 1, by minimizing the following error: finding a = 0.2 and b = 0.08, that completely defines the y W ingredient. In Figure 6, we report a comparison between the exact KE enhancement factor, the VW one and the y W , for Ar and Rn atoms. The ingredient y W is accurate near the nucleus, while the second-order gradient expansion (GE2) enhancement factor (F GE2 s = 1 + 5s 2 /27 + 20q/9) diverges at the nucleus (as q → −∞). Similar good results have been found for other noble gas atoms (data not reported). Note that y W is not a good model for r/R > 1.2, where instead, GE2 is closer to the exact kinetic enhancement factor.
The non-local ingredient y W opens the possibility of constructing KE approximations of the form: and the functional derivative of such an expression can be computed starting from its definition. After some algebra, we obtain: ∂τ ∂y W + d 3 r ∂τ ∂y W (r ) a ρ 2/3 (r ) e −bτ W (r)|r−r | 5 f (r, r ) (51) with: f (r, r ) = b 2 |∇ρ(r)| 4 |r − r | 10 32ρ 2 (r) δT W δρ(r) − Equation (51) differs from a semilocal expression mainly due to an extra integral, but it is still numerically feasible, having the same computational cost as the total KE T s .

B. Exchange Energy
For the exchange case, we start again from the observation that the 1SM model of Equation (43) is almost exact; see Figure 3. Then, instead of considering the Hartree potential of the density of the 1s-shell, we consider a screened Hartree potential of the total density. Thus, in order to describe XE density in a region near the nucleus, we introduce a β-family of ingredients, which is invariant under the uniform density scaling, of the form: where a, b and β are positive constant. Note that the exponential dumping factor of Equation (53) uses the Fermi wavelength of Equation (46). Near the nucleus of a many-electron atom, α is small, but nonzero, being a measure of the p-electrons contribution.
We recall that for any one-and two-electron systems, α = 0, the 1SM model of Equation (43) is exact, and thus, F exact x = x u . Clearly, an exact F x for one-and twoelectron systems is a very important condition in DFT that cannot be reached by any semilocal functional. In this respect, the here proposed ingredient x u (r) can be useful in further developments of non-local functionals.
In Equation (53), for simplicity, we chose β = 1, and the parameters a and b have been optimized for the Ne atom, using the same procedure as in the kinetic case. We find that for a = 0.8, b = 0.85, x u (r) is suitable for the description of the nuclear region up to r/R ≈ 0.6 − 0.8, as shown in Figure 7, in the cases of Ne and Kr atoms.
As in the KE case, non-local XE functionals with improved nuclear behavior can be developed, having a general expression (for spin-unpolarized systems): and using the generalized Kohn-Sham method 81 , the XE potential can be computed in the same manner, as has been shown in Equation (51), and will differ from a regular meta-GGA potential only due to an extra integration. Finally, we remark that the exact nuclear behavior can be reproduced only by high level functionals, which include the exact-XE density, such as the optimized-effective potential (OEP) method [82][83][84] or hyper-GGAs methods [85][86][87] , which are significantly more expensive than the proposed Equation (54).

IV. CONCLUSIONS
In conclusion, we have investigated the behaviors of the kinetic and exchange energies near the nucleus region. By employing the simple, but very powerful hydrogenic orbital model system, we have reported the exact expression for the kinetic and exchange enhancement factor near the nucleus, from the helium isoelectronic series to the semiclassical limit of a neutral atom with an infinite number of electrons. This analytical study has also proven that the 1s-model is very accurate for the kinetic (due to a subtle error compensation mainly between the 2s and 2p electrons) and for the exchange energy density.
The physics of the kinetic and exchange energy densities near the nucleus region has fully non-local features, and thus, it cannot be captured by the usual semilocal ingredients. For this reason, semilocal exchange-correlation and kinetic functional approximations are not accurate in this region. We propose density-dependent ingredients (i.e., y W of Equation (48) and x u of Equation (53)) that can well describe this important density region and can become attractive tools for the future development of DFT functionals.