The local kinetic energy density revisited

Kinetic energy density (KED) plays a fundamental role in density functional theory, chemical bonding, etc. As a basic quantity, its definition is ambiguous, which will inevitably exert unpredictable adverse effects on application. We derive another form of local total energy density (ED) and KED based on the perturbation theory. Contrary to popular belief, our study reveals that energy is only distributed inside the potential well and exactly fills the entire potential well. It seems that the ambiguity of the quantum definition of total ED is eliminable by the stability of energy distribution, and there is a one-to-one correspondence between this form of KED and total ED. As a result, the ambiguity of the KED is also eliminable. Moreover, this new form of local KED ‘selectively’ agrees well with the key regions of other commonly used ones. However, it is much less localized than other forms. Its locality features imply that it may provide a native way to describe the strong correlation and the van der Waals interactions.


Introduction
Kinetic energy density (KED), directly or indirectly, plays a fundamental role in density functional theory (DFT), chemical bonding, etc. Its role in chemical bonding has long been recognized [1], and further confirmed recently [2,3]. In DFT, kinetic energy is the leading source of non-locality [4], and the inclusion of the KED makes the M06L meta-generalized gradient approximation [5] suitable for van der Waals interactions [6,7]. For the strongly correlated systems, the potential energy dominates over kinetic energy, and the Kohn-Sham DFT that the kinetic energy functional is minimized alone to carry out the minimization of total energy functional does not work well. Giorgi et al adopted an opposite strategy, that the electron-electron repulsive potential energy is minimized alone, to handle strongly correlation and received some success [8]. The KED is basic and important, but its quantum definition has been ambiguous all the time [9][10][11]. Its definition can be roughly classified into local, semi-local and non-local. Here we focus on the local KEDs. In 1979, Cohen summarized several local KEDs that have been frequently used in literature [12]. By applying the phase space formulation of quantum mechanics, he provide us with an infinite number of local KEDs. The research and development of KED are inseparable from its application in chemical bonding [11,13,14]. In the literature, different authors prefer different forms of KED, and in general, there is no consensus on which one has absolute advantages [9]. In 2010, Anderson et al analyzed the ambiguity of local KEDs [9]. Their assessment of 'how ambiguous the KED really is?' is 'very ambiguous'. The nature of this ambiguity remains elusive. This ambiguity also extends to orbital-free DFT [15], subsystem DFT [16,17], etc. In DFT, the exchange-correlation density functional and KED functional are both unknown [18]. At present, there are some good approximation schemes for the exchange-correlation density functional [18,19]. People have also put massive efforts into the approximation functional of KED, but found it is more difficult [18,19]. In this paper, we present another form of local energy density (ED) and KED, together with some comparisons with other commonly used forms. Atomic units (a.u.), with = m = e 2 = 1, is adopted throughout this work.

Methods
Consider a 1D non-degenerate stationary perturbation problem. When there is no perturbation, it has where H 0 is the Hamiltonian of arbitrary 1D bound system with an assumption that the potential field is uniformly segmented, namely it is a piecewise constant potential [20]. We suppose the perturbation is defined by where v c is a constant, x i is the position of the ith segment, and Δx is the step size of the segmented potential field. If only take the first-order approximation of the perturbation theory, the amount of energy change can be approximately written as Now we introduce a way to strictly solve the energy change ΔE. Assuming v c is negative, then the introduce of perturbation H would form a new small potential well based on the original unperturbed potential well with the depth of v c and the width of Δx. Now, we suppose that this new small potential well is not formed by once, but in m steps. It is going to reduce a depth of δh = v c /m from the previous potential at each step, then v c is divided into m tiny δv c layers. For each step, δv c can be regarded as a perturbation of the previous potential, and the perturbation method can be applied continuously to this process, then where Ψ j n is the wave function at the jth step, namely the wave function after introducing the superimposed perturbations jΔh/m. Noting that all of the potential wells are piecewise uniform, so the corresponding Schrödinger equation of each step can be nearly solved exactly with the transfer matrix method [20,21]. Let Δh → 0, then The approximately equals sign in formula (5) is eliminated, and v c is no longer required to be a perturbation. The energy change of the system is given strictly by formula (6). Since its integral is only performed inside the H potential well, it implies that the ΔE part of energy maybe only distributes inside the H potential well. Then, the integrand |Ψ n (h, x)| 2 is relevant to the ED. The square of wave function represents the probability density, but here |Ψ n (h, x)| 2 is proportional to the energy area density, just dimensionally different. Let H be the dimensions of h, namely, the dimensions of energy. Then |Ψ n (h, x)| 2 H should denote energy area density. For convenience, we will directly call |Ψ n (h, x)| 2 the energy area density in the following. The above results can be generalized to the arbitrary potential well V(x). We could let the original potential field to be a zero field, and take H = V(x). Perform the above layering operation on the potential well V(x), then The upper and lower integral limits h b and h 0 are the potential energy at well bottom and well mouth, respectively. |Ψ n (h, x)| 2 is the energy area density of the single nth orbit. The total energy of state Ψ n is We will set the well mouth as potential zero in the following numerical calculation. The energy line density of single orbit at position space is given by When 2N noninteracting fermions occupy N orbits, then In the above derivation, it is assumed that the potential well V(x) is stepwise formed. That is, V(x) is 'horizontally' divided into m layers, then the potential well at the jth step is where h j = jΔh/m. Actually, it is also implicitly supposed that the subsequent layering does not affect the previous energy distribution. If the energy distribution remains unchanged during stepwise layering, the potential energy should do the same. Therefore, the potential energy for a specific energy level in the jth where E j and E j−1 are the total energy at the jth and (j − 1)th step respectively. The distribution function of δE j part of energy is |Ψ j (x)| 2 , which is the probability itself. Similarly, the distribution of potential energy in the jth layer should be proportional to probability, that is δV j |Ψ j (x)| 2 /δE j . So, the distribution of kinetic energy in the jth layer is |Ψ j (x)| 2 − δV j |Ψ j (x)| 2 /δE j . Then the KED of the entire potential well is Specifically, Ψ j is the eigenvector of the Schrödinger equation with potential energy to be V j in formula (11), and δV j , δE j are the potential energy and the eigen energy in the single jth layer respectively. δh j is the thickness of the jth layer. As claimed by the constraint condition of summation, the distribution is restricted to the inside of the potential well, namely only the contribution within layers is considered. For convenience, the energy level index n in equation (12) has been omitted. Since δV j < δE j is always true, so K F is positive definite. The entropy calculated by K F is also positive definite [22,23]. The average kinetic energy K F satisfies the Virial theorem and is satisfied within each layer. Please refer to the appendix for more detailed derivation, numerical verification and parameters for the following numerical results. We chose several commonly used forms of local KED to compare with K F . As summarized by Cohen [12], they are and K D can also be expressed as −( 2 /2m)( 1 4 ∇ 2 |ψ| 2 − |∇ψ| 2 ). In this work, the above ψ is the same as Ψ j=m . K A and K B are also called the Schrödinger form and Kohn-Sham form of KED respectively. In the literature, the potential ED is generally expressed as The corresponding total ED can be written as Figure 1 exhibits the energy area density of a uniformly segmented 1D harmonic oscillator. The jagged outline of the ED distribution is exactly the curve of the potential field. For the corresponding results of 1D hydrogen atom and finite square potential well, please refer to the end of appendix . All the energy distributes inside the potential well. It shows that the energy of the ground state mainly distributes in the deep region of the potential well, while the energy distribution of excited state is lifted. Figures 2(a) and (c) compares several commonly used forms of ED with ξ F . In the near-nuclear region of figure 2(a), ξ F , ξ B , and ξ D rise rapidly in a similar trend, but only ξ F eventually approaches zero at the   (9) and (17) respectively. K A , K B and K D are given by equations (13)-(15) respectively.

Calculation results
nuclear center, which is consistent with ξ A . Besides, the distribution range of ξ F is appreciably different from others. Figures 2(b) and (d) list the kinetic energy distributions. Firstly, the integrals of different forms of KED are consistent. In figure 2(b), K F curve has a shape that is comparable to K A , but K F has no negative density. Furthermore, K F and K B are basically coincident in the interval [2,4]. In figure 2 The degree of locality of K F significantly differs from others. Particularly, in the long-range area of Coulomb potential, K F and ξ F decay much more slowly than other forms, as shown in the inset of figures 2(a) and (b). This feature should be beneficial to describe the long-range interaction between layers. By analyzing the K F distribution for other more types of potential fields, we obtain the general rule of the locality. In deep wells (strongly bound), K F distributes wider than others; in shallow wells (weakly bound), K F is more localized relatively. For square well, the degree of locality of K F and K D is always the same.
The commonly used forms of KED are not universal, but from a partial perspective, it is conceivable that some KED forms can reproduce the kinetic energy distribution in a certain region of a specific scenario, such as near-nuclear and middle-range regions. As shown in figure 2(b), the trend of K B and K D near the nucleus is opposite to that of K A and K F . Chattaraj and co-workers found evidence in the entropy contours of water molecules that K B does not reproduce the KED near nucleus [23]. So, K A and K F ought to be more representative of the real distribution in this scenario. The K B form is widely used in chemical . The KED and ED for multiple orbits occupied by noninteracting fermions in harmonic oscillator potential. N is the number of orbits that filled by 2N fermions. ξ F is given by formula (10). ξ A,B,D is given by formula (17), with [26]. To accommodate more orbits, the potential well is truncated at x = ±7.5 a.u. in numerical calculation.
bonding [10,[27][28][29][30]. The bonding interval is mainly in the medium-range (∼1-4 a.u.) [31]. As the KED of the Coulomb potential shown in the inset of figure 2(b), K B and K F are basically coincident in range 2-4 a.u. In addition, as for the harmonic oscillator in figure 2(d), the shape of K F and K D are quite similar. Accordingly, it is possible that the forms of KED that used in various scenarios can be summarized into K F . Figure 3 shows the ED and KED distribution curves of noninteracting fermions in 1D harmonic oscillator. When N is large enough, different forms of the total ED tend to be consistent, so is the KED. The K D curve is always smooth and therefore has received special attention [26]. But we notice that the corresponding total ED ξ D fluctuates more markedly than ξ A . In contrast to classical particles, we incline to believe that the fluctuation of total ED should not be greater than that of KED. At N = 8, ξ F is the smoothest one, with large fluctuations only around its maximum. When N = 25, ξ F exhibits large fluctuation in the range [−2, 2]. Its amplitude is comparable to ξ A , and is significantly smaller than that of ξ B and ξ D . The fluctuation of the corresponding K F is stronger than the other forms. There is a visible difference between K F and others. K F first rises extremely rapidly near the edge of the potential well, and then shows wave properties, as shown in figure 3(a). The behavior of kinetic energy rising rapidly at the edge of the harmonic oscillator potential is identical to that of the classical particles.

Discussion
Equations (6) and (7) can be derived from the perturbation theory. However, a prerequisite must be satisfied if |Ψ(h, x)| 2 is to have any real physical meaning. That is, when v c or V(x) is stepwise layered, the energy distribution of the previous potential field remains unchanged. In fact, there are infinite schemes for dividing (or meshing) the potential field, and each scheme corresponds to a form of energy distribution, but not all of the schemes can meet the above prerequisite. We choose the 'horizontal' layering scheme. We argue that this layering scheme can keep the energy at the same depth stable, thus maintaining the overall stability. The ambiguity in total ED is eliminated by the stability of energy distribution. The distribution of the total ED is proportional to the probability. Meanwhile, for a certain layer, the potential field is uniform, so it is natural to assume that the potential ED inside the layer is also proportional to the probability. Based on it, we derived the KED, which is again proportional to the probability of each layer. Similarly, it was proved by Nagy that if let the local temperature be a constant for the whole system, the KED would be proportional to the electron density, and the information entropy would be maximized [32]. On the whole, it should be the most reasonable and natural to conclude that the potential ED inside the layer is proportional to the probability. As there is a one-to-one correspondence between this form of KED, potential ED and total ED, it can be further considered that the ambiguity of KED is also eliminable.
Whether the above distribution is true or not is closely tied with a simple question. After introducing the perturbation H (Assuming that H locates at the bottom of the original potential well), the total energy change of the system is exactly given by the expression (6). Then, does this part of energy distribute in the Δx interval where H is located? Or does it distribute throughout the position space? Both cases can be explained accordingly. The former is a direct interpretation of formula (6). The latter obeys the conventional view.
In addition, this form of distribution matches a simple physical picture for the perturbation theory, and it can be used to simplify the perturbation calculation, as detailed in the appendix .
The above results can be extrapolated to the cases of three-dimension and many-body. The definition of KED we derived is actually the superposition of KEDs of a series of quantum states. From the perspective of interaction, these quantum states include intermediate states from non-interacting states to actual interacting states. The interaction of these intermediate states is weakened, and the degree of locality of the KED of these weakened states is correspondingly reduced. Then the locality features should be consistent with the 1D situation at the cases of higher dimension and many-body.
Due to its locality features, K F is potentially valuable for dealing with strong correlation and van der Waals interactions. Briefly, the failure of traditional Kohn-Sham DFT in these two cases can be attributed to that the interaction distance are long relatively [8,33,34]. It has been shown that the KED is responsible for interactions [2,3,6,10]. However, the Kohn-Sham KED K B and others are too localized to describe long-range interactions. To conquer the difficulties, people had directly or indirectly introduced the semi-local and non-local functionals [34][35][36][37]. K F may meet the non-local requirements for describing long-range interactions, and consequently has a good application prospect in DFT.

Summary
Here another pattern of energy distribution is derived. Unconventionally, all of the energy is distributed inside the potential well and exactly fills the entire potential well. Some comparisons with other commonly used forms is made. The KED K F agrees well with the key regions of other forms. But, the locality of K F significantly differs from others. This makes it have a fantastic prospect for describing long-range interactions in a native way. In our method, the perturbation itself is regarded as a differential element. Through integration, this method can be used to deal with any non-perturbation problem with the perturbation theory. This work is inspired by the result of a numerical experiment. Here we will start with this numerical experiment. When numerically solving the Schrödinger equation, it is often necessary to discretize the potential field. Generally, the discretized potential is an approximation of the original potential and introduces errors. The potential error at different regions has different effect on the accuracy of the solution. For the same size error, the error far away from the center of the potential field has less influence on the solution. The mesh of this area does not need to be too dense, such as a dense mesh is usually not required near the boundary. The error in the central region of potential field has a greater impact on the accuracy of the solution, and a more dense mesh is needed.
To study the influence of potential error at different regions on the accuracy of the solution, we designed a numerical experiment. For simplicity, we consider the 1D potential well problem. Consider a piecewise uniform, namely piecewise constant, potential V(x), such as is a smooth function, and x i+1 = x i + Δx, as shown schematically in figure A1. Since V(x) is piecewise uniform, the transfer matrix method [20,21] can solve the bound states of the potential V(x) with high precision. To investigate the influence of potential error at different spatial positions, we adjust the potential field at different positions successively. When taking turns to the spatial position x i , the current potential field becomes (A.1) Figure A1. The traditional 1D harmonic oscillator potential is discretized to a piecewise potential, and truncated at x = 1.5, with the step size Δx = 0.5. The vertical axis is the numerical value of the curves in atomic units. V(x) and V (x) are the potential without and with perturbation respectively. In this case, In other words, the perturbation forms a new potential well based on the original potential well, as the square area colored in orange. The arrow marks the ED distribution at the upper and lower edges of H . For the ground state, the energy change caused by H , or the energy filled in v c well, can be approximately expressed as: , where E and Ψ are the eigenenergy and eigenvector for v c = 0 respectively. Similarly, E i is the normalization of the influence of error on eigenenergy. The curves of E i − |Ψ| 2 reflects the relationship between E i and probability density |Ψ| 2 . When v c is very small, they are almost equal.
v c is a constant, and can be regarded as a man-made error which can also be viewed as a perturbation. Since v c is constant, V i (x) is still piecewise uniform, whose eigenstate can also be solved precisely. Thus, the influence of v c on the solution can be accurately observed.
Here, we only consider the ground state. The eigenenergy E i of potential field V i (x) with different i are solved successively, and the E i (x i ) curves are plotted in figure A2, which directly reflects the influence of potential error at different regions on eigenenergy. We find an interesting phenomenon in figure A2(a) that the E i (x i ) curves have considerable similarity to the probability density |Ψ(x i )| 2 . The influence of the error v c on the eigenenergy is nearly proportional to the probability density. When v c is tiny, its effect on eigenenergy is small correspondingly, but the linear transformation E i (x i ) almost coincides with probability density |Ψ(x i )| 2 as shown in figure A2(b), indicating that there may be some equality relationship between the probability density and the change of eigenenergy. Figure A3(a) shows the relationship between the energy change ΔE and v c , where ΔE = E i − E, and E is the exact eigenenergy of the ground state when v c = 0. Within the range shown, ΔE(v c ) of different positions are all linear. We find that the slope value of ΔE(v c ) at different spatial positions is proportional to its probability density, which means that ΔE should be proportional to v c , and the proportionality coefficient is related to probability. Since the magnitude of ΔE is 10 −3 , which is a small amount, the man-made potential error v c can be regarded as a perturbation. If taking the first-order approximation, Since v c is a perturbation, its influence on the wave function is very limited. As a result, the probability amplitude can be approximated as unchanged after the perturbation is added. So when v c is small, ΔE is approximately a linear function of v c , which has explained the linear relationship between ΔE and v c in figure A3(a). Since the integrating interval is [x i , x i+1 ], there are reasons to assume that the energy represented by ΔE fills in the region where v c is located. This integral expression of energy could relate to the ED [24,38,39], and v c |Ψ| 2 could be connected with the energy distribution. The expression x i+1 x i v c |Ψ| 2 dx is an approximation of ΔE, and it is only suitable for the case of perturbation. Now we introduce a way to strictly solve the amount of energy change ΔE, regardless of whether v c is perturbation or not. Assuming v c is negative, then the introduction of v c would form a new potential well based on the original potential well V(x) with the depth of v c and the width of Δx. Now, we suppose that this new potential well is not formed by once, but in m steps. It is going to reduce a depth of δh = v c /m from the previous potential at each step, as shown in figure A4, then v c is divided into m small δv layers. For each step, δv can be regarded as a perturbation of the previous potential, and the perturbation method can be applied continuously to this process, then Ψ 0 is the wave function when δh is not introduced, Ψ j is the wave function at the jth step with the current potential V i (x) expressed as .
(A.5) When δh → 0, namely m → ∞, then The integral interval of expression (A.6) is limited to the inside of the potential well. |Ψ(h, x)| 2 is proportional to the area density of energy, just dimensionally different. For convenience, we directly call it as the energy area density. Now, with the view of energy area density, we go back to the relationship of ΔE and v c . If the depth of the v c well is small enough, and with the above analysis that there is energy area density inside the v c well, its ED along the depth direction is approximately uniform. Let E p be the energy that fills in the potential well v c . Then it can be approximated as which is consistent with the result of the perturbation theory. As shown in figure A1, this may be a physical picture for the perturbation theory, which can be used to optimize perturbation calculation. Follow the linear approximation detailed in the next paragraph, the amount of energy change after the introduction of perturbation H can be more precisely expressed as where Ψ (0) and Ψ (1) are the zero-order and the first-order approximation of wave function in perturbation theory respectively. Ψ (1) can be replaced with a higher approximation for better accuracy, as shown in figure A3(b). It shows that the linear approximation is better than the second order approximation of perturbation theory. The variation of ED in the depth direction should be continuous. For equation (A.4), it has approximatively assumed that the filling of the ED inside each δv single layer is uniform. Undoubtedly the Figure A5. Generalize the layering scheme to the whole potential field. Start with a zero potential field, the corresponding wave function of the ground state equals to zero, which means that the ED distribution at the well mouth is zero. h 1 , h 2 and h 3 in panel (a) are the heights of the first, second and third stair respectively. In this case, the total energy of the ground state can be approximately written as E ≈ The more the number of layers, the more accurate the result. linear approximation would be better. As shown in figure A4, the energy of each layer can be written as H Jiang Figure A7. The ED and KED for single particle at the ground state corresponding to the potential type of figure A6. and the change of total energy ΔE can be approximated as where the expression of energy change with linear approximation is labeled as E L . Three different approximations of energy change influenced by v c are plotted in figure A3(b), which shows the relationship between the energy change for different approximates and the exact energy change. The linear approximation E L is generally better than the first-order and second-order approximations of the perturbation theory. The first-order and the second-order approximations of λ in the perturbation theory are both the first-order approximation of v c . Only in a few cases, the second-order approximation of perturbation theory is the second-order of v c , such as the interval [−10, −0.1] in the upper right panel of figure A3(b). In the cases listed, the linear approximations E L are all second-order of v c . When v c is small, all value of ΔE/E p for different approximations have a tendency of approaching 1, which indicates that the amount of energy filled in v c is exactly the decrease of system energy. When v c is large, the error of the normal perturbation approximation is large. But, the layering method described above is suitable for any size of v c . Figure A3(c) plots the relationship of E L and ΔE. Even if v c is no longer a perturbation, E p calculated by expression (A.10) can be sufficiently close to the exact ΔE. It again confirms the expressions we derived, namely the change of eigenenergy influenced by v c can be exactly calculated by expression (A.6).
The width of the v c well discussed above is just right the step size Δx. When the width of the well is arbitrary, similar stepwise accumulation calculation can also be performed. To generalize it to the arbitrary potential field V(x), we can suppose that the potential field V(x) is formed from a zero field through m steps, and the potential field is reduced by a part of depth at each step, as shown in figure A5. Calculate the wave function at each step and integrate them inside the potential well, then the expression of ΔE is When m is large enough, E L approaches to the exact total energy E. Then the expression (A.11), namely the expression (7) in the manuscript, is verified. Our derivation is based on the perturbation theory, but it differs from the conventional perturbation theory, mainly in two aspects: 1) Here the wave function with the perturbation is exactly knowable, so the wave function we used does not have the problem of convergence. 2) We take the limit of the perturbation, and integrate it. In the case of taking the limit, the high-order approximation can be ignored. We finally get the exact expressions. They are applicable in the case of introducing perturbation at the zero filed. From the view point of energy area density, the ED distributed at the border of the layer is always accurate, while the energy distribution inside the layer is linearly approximated. In this way, even if the thickness of the layer is