#  A general and systematic approach for kinetic energy density functional:
Green’s function formalism.

Priya  Department of Chemistry, Indian Institute of Technology Kanpur, Uttar
Pradesh, India  Mainak Sadhukhan<EMAIL_ADDRESS>Department of Chemistry,
Indian Institute of Technology Kanpur, Uttar Pradesh, India

###### Abstract

An accurate expression of the kinetic energy density of an electronic
distribution in terms of a single particle reduced density matrix for
atomic/molecular systems is a long-standing problem in electronic structure
theory. Existing kinetic energy density functionals are based on model
potentials and require empirical parametrizations to fit the accurate values
of the kinetic energy of atoms. Moreover, the correct kinetic energy density
which produces proper atomic shell structure is still an unsolved problem. In
this work, we have developed an exact methodology which can be used to
systematically derive kinetic energy of an electronic system. We have also
shown how the kinetic energy is dependent upon number of electrons $N$ and the
atomic charge $Z$ using a qualitatively similar model system. We have also
connected our methodology with other forms of Green’s functions.

##  I Introduction

An accurate approximation of kinetic energy density (KED) functional is the
main bottleneck in the development of orbital-free density functional theory
(OF-DFT)Ludeña and Karasiev ; Wesolowski and Wang (2013); García-Aldea and
Alvarellos (2008); Perdew and Constantin (2007); García-Cervera2008968;
Ghiringhelli and Delle Site (2008); Eek and Nordholm (2006); Karasiev  _et
al._ (2006, 2009); Trickey  _et al._ (2009); Huang and Carter (2010a);
Karasiev and Trickey (2012). In this work, we have developed a general, spin-
multiplicity-independent formalism to systematically calculate the KED for any
Fermionic system in its ground electronic state. Our formalism does not
recourse to any model-based KED as the zeroth order approximation.

One particle reduced density matrix $\rho(\mathbf{r})$ is the fundamental
variable in DFT. Since $\rho(\mathbf{r})$ is a three-dimensional quantity
irrespective of system size, it should be obtained from a single Euler-
Lagrange equationHohenberg and Kohn (1964)

$\frac{\delta E[\rho]}{\delta\rho}=\mu$ (1)

where $E[\rho]$ and $\mu$ are total electronic energy as a functional of
$\rho(\mathbf{r})$ and chemical potential respectively. The functional form of
total energy $E[\rho]$ is unknown rendering the Eq.(1) exactly unsolvable. The
$E[\rho]$ contains potential energy $V[\rho]$ and kinetic energy $T[\rho]$. In
contrast to $V[\rho]$, for which adequately accurate approximations exist,
similar approximations are yet to be attained for $T[\rho]$. Note that
$V[\rho]$ can be broken down into classical Coulomb interaction $J[\rho]$ and
non-classical exchange-correlation (XC) energy $E_{XC}[\rho]$. However, no
such universally agreed-upon breakdown of $T[\rho]$ exists in general. Apart
from an accurate kinetic energy value, $\frac{\delta T[\rho]}{\delta\rho}$
should give characteristic features of the electronic system such as shell
structure for atomic and molecular electron density and the dependence of KED
on number of electrons $N$ and atomic charge $Z$. Kohn-Sham DFT (KS-DFT) Kohn
and Sham (1965) bypasses this challenge by introducing $N$ orbitals and
merging the resulting kinetic energy corrections in exchange-correlation
functionals. While KS-DFT is in principle an exact approach, the computational
cost increases significantly with the number of particles due to the need to
solve $N$-coupled differential equations. As a result, the quest for OF-DFT
and hence an accurate form of $T[\rho]$ remains a challenging yet essential
problem.

The problem of finding an accurate approximation of KED is related to its non-
local definition in terms of the one-particle density matrix (ODM)
$\rho_{1}(\mathbf{r},\mathbf{r}^{\prime})$. The kinetic energy and the KED
$t(\mathbf{r},\rho(\mathbf{r}))$ of an electronic system are defined by
(atomic units have been used throughout the article unless mentioned
explicitly)

$T[\rho]=\int t(\mathbf{r},\rho(\mathbf{r}))\text{d}\mathbf{r}$ (2)

and

$t(\mathbf{r},\rho(\mathbf{r}))=-\frac{1}{2}\nabla^{2}\rho_{1}(\mathbf{r},\mathbf{r}^{\prime})|_{\mathbf{r}=\mathbf{r}^{\prime}}.$
(3)

All existing methods to compute the KED start with approximations of
$\rho_{1}(\mathbf{r},\mathbf{r}^{\prime})$Parr and Yang (1994).

The earliest approximation of KED was modeled after homogeneous electron gas
(HEG) Thomas (1927); Fermi (1927) followed by von Weizsäcker’s correction for
slightly inhomogeneous electron density Weizsäcker (1935). The sum of HEG
kinetic energy (Thomas-Fermi KE)

$T_{TF}[\rho]=C_{TF}\int\rho^{5/3}(\mathbf{r})d\mathbf{r}$ (4)

and von Weizsäcker KE

$T_{vW}[\rho]=\frac{1}{8}\int\frac{\left(\nabla\rho(\mathbf{r})\right)^{2}}{\rho(\mathbf{r})}d\mathbf{r},$
(5)

however, overestimated the kinetic energy of the system. Note that the
$T_{vW}$ is the exact kinetic energy for electrons in $1s$ orbital of a
hydrogen atom. Moreover, they are inadequate to describe molecular binding
Teller (1962). Due to this overestimation a scaled von Weisäcker term in the
form

$T[\rho]=T_{TF}[\rho]+\lambda T_{vW}[\rho]$ (6)

has been considered, known as generalised Thomas-Fermi-von Weizsäcker kinetic
energy (TF$\lambda$W)S. Clearly, $\lambda=1$ represents the upper limit of KE.
It has been found empirically that $\lambda=1/5$ is the best fit for Hartree-
Fock electron density. Analytical derivation of $\lambda$ using generalized
gradient approximations (GGA) in turn predicted $\lambda=1/9$. This
discrepancy has not been resolved and no shell structure is obtained in
electron density using TF$\lambda$W theory. As a consequence, standalone
quantum chemical techniques could not be developed based on this approach Wang
and Teter (1992a). Using the $T_{TF}$ as the zeroth order approximation is the
suspect here for the realistic atoms are more closely represented by a
hydrogen atom than HEG.

In the condensed matter physics and materials science community, the kinetic
energy density functional (KEDF) is enhanced using non-local functionals to
improve upon the Thomas-Fermi (TF) and von Weizsäcker (vW) terms. The Wang-
Govind-Carter (WGC) functional is a non-local KEDF that incorporates density-
dependent terms and gradient corrections and is defined as
$T_{WGC}[\rho]=T_{TF}[\rho]+T_{vW}[\rho]+T_{NL}[\rho]$. Witt  _et al._
(2018a); Wang  _et al._ (1999); Carter_Lignères2005; Hung and Carter (2009).
Non-local functionals, such as generalized gradient approximation (GGA)-type
functionals and density-dependent kernels, are also employed to improve the
local TF model Wang and Teter (1992b); Huang and Carter (2010b); Mi  _et al._
(2018); Shao  _et al._ (2021); Carter_Lignères2005; Witt  _et al._ (2018b).
Recently, significant work has been done using machine learning to find KEDFs
for solids, atomic, and molecular systems Snyder  _et al._ (2012); Meyer  _et
al._ (2020); Sun and Chen (2024); Seino  _et al._ (2018); Hollingsworth  _et
al._ (2018). Another important term related to KED is the Pauli potential
($T_{P}$) Tal and Bader (1978); Bartolotti and Acharya (1982); Levy and Ou-
Yang (1988), which is the difference between the KS kinetic energy ($T_{s}$)
and the von Weizsäcker kinetic energy ($T_{vW}$), defined as
$T_{P}=T_{s}-T_{vW}$. The Pauli potential is often approximated using
bifunctional expressions and is used in approximating molecular Pauli
potentials Finzel (2021a) as well as in quantum crystallography Finzel
(2021b). Additionally, the Pauli potential is approximated for spin-polarized
atoms.

In contrast, for closed-shell systems, described by a single determinantal
wavefunction the kinetic energy can be shown to be described by the sum of
full von Weizsäcker term and a correction term Alonso and Girifalco (1978).
This correction term can be computed from a correlation function. It has been
shown from this formalism that the KE and XC energy are related through joint
density matrices, and therefore are related to the aforementioned correlation
functionAlonso and Girifalco (1978); Alonso and Balbás (1981); Deb and Ghosh
(1983).

The traditional kinetic energy functionals does not give the characteristics
features of electron density like the shell structure in atoms. A class of
work based on quantum-fluid dynamical density-functional approach Deb and
Ghosh (1982), introduces two functions $f(\mathbf{r})$ and $g(\mathbf{r})$ for
an accurate representation of electron densities in many-electron system where
they are modeled after Hartree-Fock densitiesDeb and Ghosh (1983). This
ensures not only the accuracy of kinetic energy but also gives the shell
structure. While promising, this line of approach is limited to closed shell
systems and requires a model-based form of correlation function. Deb-Chattaraj
equation is an hydrodynamical equation used to model the dynamic electron
densities in many-electron systems Deb and Chattaraj (1989).The Time-dependent
hydrodynamical equation (TDHE) is derived from the Schrödinger equation within
a quantum hydrodynamical framework. This approach enables the direct
computation of ground state electronic densities, eliminating the need for
iterative self-consistent procedures (SCF). It has been applied to calculate
ground state properties for noble gas atomsRoy  _et al._ (1999).

Non-local approximations to correlation function have been considered
subsequently Alonso and Balbás (1981).The correlational functional is defined
as
$C(\mathbf{r},\mathbf{r}^{\prime})=\frac{\Gamma(\mathbf{r},\mathbf{r}^{\prime})}{\rho(\mathbf{r})\rho(\mathbf{r}^{\prime})}-1$.

$T[\rho]=T_{vW}[\rho]+C_{TF}\int
f(\mathbf{r})\rho^{5/3}(\mathbf{r})d\mathbf{r}$ (7)

for bare $T_{TF}[\rho]$ the $f(\mathbf{r})=1\forall\mathbf{r}$. Had the
electron been bosons, all of them would occupy the ground state $1s$ orbital.
In such situation, total kinetic energy will be completely described by
$T_{\text{vW}}[\rho]$ implying $f(\mathbf{r})=0\forall\mathbf{r}$. Therefore
$f(\mathbf{r})$ enforces the Pauli exclusion principle as well as the effect
of inhomogeneity in kinetic energy. Since the shell structure is a direct
consequence of the Pauli exclusion principle, $f(\mathbf{r})$ manifests the
fermionic nature of electrons in kinetic energy. As discussed above the
$f(\mathbf{r})=1$ is the Thomas-Fermi term which is similar to
$\hat{\mathcal{H}}_{Z}\to 0$, in atoms away from the nucleus as the electron-
nuclear interaction part tends to zero i.e. and electron-electron repulsion
becomes dominating interaction.

Apart from being specific model-dependent these previous works also fail to
explain how $T[\rho]$ and $t(\mathbf{r},\rho(\mathbf{r}))$ are dependent on
$N$, $Z$, and the spin multiplicity of the system. Also, the accuracy of the
approximation involved was bench-marked against the global (integrated) value
of $T[\rho]$, which fudges local details of the
$t(\mathbf{r},\rho(\mathbf{r}))$ which is quintessential to produce correct
shell structure. Therefore, a systematic, specific model-independent
methodology is required which should be able to explain $N$ and $Z$ dependency
and can be employed for arbitrary spin-multiplicity.

In our recent work Priya and Sadhukhan (2022), we developed a systematic
methodology to find kinetic energy density, using Green’s function (GF) which
satisfies all aforementioned qualities. We will call it Green’s function
formalism (GFF). The GF employed in this present work is defined as the
Laplace transform of the ODM

$G(\mathbf{r},\mathbf{r}^{\prime};\beta)=\int_{0}^{\infty}\text{d}\epsilon_{F}e^{-\beta\epsilon_{F}}\rho_{1}(\mathbf{r},\mathbf{r}^{\prime})\equiv\bra{\mathbf{r}}e^{-\beta\hat{\mathcal{H}}}\ket{\mathbf{r}^{\prime}}.$
(8)

Here $\epsilon_{F}$ and $\hat{\mathcal{H}}$ are Fermi energy and Hamiltonian
of the system respectively. The auxiliary variable $\beta$ is mathematically
similar to the inverse temperature used in Statistical mechanics. However, in
this formalism the similarity is only notational and does not signify any
physical thermodynamic temperature. The present formalism is a $0$ K theory.
The ODM can be retrieved from the GF via inverse Laplace transform Parr and
Yang (1994).

$\rho_{1}(\mathbf{r},\mathbf{r}^{\prime})=\lim_{\Gamma\to\infty}\frac{1}{2\pi
i}\int_{\gamma-i\Gamma}^{\gamma+i\Gamma}\frac{\text{d}\beta}{\beta}e^{\beta\epsilon_{F}}G(\mathbf{r},\mathbf{r}^{\prime};\beta)$
(9)

where $\gamma\in\mathbb{R}_{>}$. In this work, we have developed a systematic
expansion technique for the GF applicable for atomic and molecular electron
densities akin to perturbation theory.

After a brief introduction of existing approaches in Section I, we formulate
our methodology in Section II.1. This exposition is followed by a discussion
regarding the connection between traditional Green’s functions (GF)Economou
(2006) and our Green’s function formalism (GFF) in Section II.2. To
demonstrate the efficacy of the GFF methodology, we applied it to simple model
systems, such as a one-dimensional harmonic oscillator in the weak interaction
limit in Section III.1 followed by a comparison with the results from
traditional GF in Section III.2. After analyzing the spatial nature of the GF
in Section III.3 we analyze the dependence of kinetic energy on the number of
electrons ($N$) and the Fermi energy ($\epsilon_{F}$) in Section III.4. After
testing the methodology on noble gas atoms (Section III.5) we show how a limit
of the weak interaction potential ($W_{0}$) naturally arises from our
formalism in Section III.6. The paper concludes with Section IV.

##  II Analytical formalism

###  II.1 Green’s function formalism

To begin our Green’s function formalism (GFF), we define the atomic
Hamiltonian

$\hat{\mathcal{H}}=\underbrace{\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{Z}}_{\hat{\mathcal{H}}_{H}}+\hat{W},$
(10)

where $\hat{\mathcal{H}}_{0}$, $\hat{\mathcal{H}}_{Z}$, and $\hat{W}$ are
free-particle Hamiltonian, electron-nucleus attraction potential and inter-
electronic repulsion potentials, respectively. The Hydrogenic Hamiltonian with
atomic charge $Z$ is defined as
$\hat{\mathcal{H}}_{H}=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{Z}$. The
corresponding GFs are

$G^{H}(\mathbf{r},\mathbf{r}^{\prime};\beta)=\bra{\mathbf{r}}e^{-\beta\hat{\mathcal{H}}_{H}}\ket{\mathbf{r}^{\prime}}$
(11)

and

$G^{0}(\mathbf{r},\mathbf{r}^{\prime};\beta)=\bra{\mathbf{r}}e^{-\beta\hat{\mathcal{H}}_{0}}\ket{\mathbf{r}^{\prime}},$
(12)

respectively. For atomic systems, the contribution of the inter-electronic
repulsion $\hat{W}$ is smaller compared to $\hat{\mathcal{H}}_{H}$ (_vide_
Eq.(10)) 111Consider Helium and Lithium atoms with exact non-relativistic
energies as $\approx-2.904$ a.u. and $\approx-14.864$ a.u., respectively. To
them the hydrogenic parts of the Hamiltonian attribute $\approx-4.0$ a.u. and
$\approx-20.25$ a.u. The contributions of the inter-electronic repulsion are
evidently smaller compared to the hydrogenic part.. As a result,
$G^{H}(\mathbf{r},\mathbf{r}^{\prime};\beta)$ is considered as the natural
leading order term of $G(\mathbf{r},\mathbf{r}^{\prime};\beta)$. This is in
contrast with previous attempts where the leading order term was chosen from
an arbitrary model system Yang (1986); Brack  _et al._ (1976). HEG Hamiltonian
was almost exclusively used as this leading order Hamiltonian despite the fact
that the atomic and molecular systems are qualitatively very different from
HEG. In Section I, we have presented the true kinetic energy of an atom as the
admixture of two extreme limits viz. HEG ($T_{TF}$) and von Weizsäcker
($T_{vW}$) kinetic energies. The $T_{vW}$ represents a hypothetical Bosonic
hydrogenic system where all particles occupy the ground state of the atom
_viz_ $1s$ orbital. Note that the hydrogenic GF $G^{H}$ employed in our
approach (Eq.(11)) is of Fermionic nature and therefore yields kinetic energy
(say $T_{H}$) which is different from the $T_{vW}$. We can relate $T_{H}$ and
$T_{vW}$ in the following manner. Consider a Bosonic system where all the
population is on the ground state. On the other hand, the Fermions follow the
Pauli exclusion principle leading to exchange-hole. The exchange-hole is
responsible for pushing the electrons to higher energy statesDyson and Lenard
(1967); Lieb and Thirring (1975). If we switch off the exchange-hole, all
electrons occupy the ground state $1s$. As a result, $T_{H}$ become $T_{vW}$.
Since the kinetic energy of the excited states is higher than that of ground
state, the latter represents a part of the former _ie_ $T_{H}>T_{vW}$.

On the same note we can interpret the interaction GF $G^{int}$ defined as the
difference between total GF $G$ and the hydrogenic GF $G^{H}$

$G^{int}(\mathbf{r},\mathbf{r}^{\prime};\beta)\equiv\bra{\mathbf{r}}\hat{O}\ket{\mathbf{r}^{\prime}}=G(\mathbf{r},\mathbf{r}^{\prime};\beta)-G^{H}(\mathbf{r},\mathbf{r}^{\prime};\beta).$
(13)

Here

$\hat{O}=e^{-\beta\hat{\mathcal{H}}}-e^{-\beta\hat{\mathcal{H}}_{H}}.$ (14)

Since $G^{int}$ contains inter-electronic repulsion, this term represents the
central challenge for our formalism.

In order to get the systematic expansion of $G^{int}$ we expand $\hat{O}$
using Zassenhaus formula Casas  _et al._ (2012); Magnus (1954) as

$\hat{O}=-e^{-\beta\hat{\mathcal{H}}_{H}}\left(1-e^{-\beta\hat{W}}e^{-\frac{\beta^{2}}{2}\commutator{\hat{\mathcal{H}}_{H}}{\hat{W}}}\ldots\right).$
(15)

We find

$\commutator{\hat{\mathcal{H}}_{H}}{\hat{W}}=\commutator{\hat{\mathcal{H}}_{0}}{\hat{W}}$
(16)

since $\commutator{\hat{\mathcal{H}}_{Z}}{\hat{W}}$ commutes. However, such
replacements cannot be done for terms containing triple commutators onwards,
leading to our first approximation

$\hat{O}\approx
e^{-\beta\hat{\mathcal{H}}_{H}}e^{\beta\hat{\mathcal{H}}_{0}}\left(e^{-\beta(\hat{\mathcal{H}}_{0}+\hat{W})}-e^{-\beta\hat{\mathcal{H}}_{0}}\right).$
(17)

Note that despite this approximation in the third-order commutator onwards,
the expansion of the second-order commutator represents a partial sum up to
the infinite order interaction. Therefore it is of the similar spirit of
random phase approximation Gell-Mann and Brueckner (1957). Additionally, this
approximation is qualitatively acceptable where
$\hat{\mathcal{H}}_{0}\approx\hat{\mathcal{H}}_{H}$, away from the nucleus
where $\hat{\mathcal{H}}_{Z}\to 0$. Note that we have not completely replaced
$\hat{\mathcal{H}}_{Z}$ by $\hat{\mathcal{H}}_{0}$ since the leading term of
the expansion is still $G^{H}(\mathbf{r},\mathbf{r}^{\prime},\beta)$. Also,
this approximation is consistent with the fact that $f(\mathbf{r})\to 1$ away
from the nucleus. We found, however, that this approximation leads to an
unphysical asymmetry in the resulting GF with respect to $\mathbf{r}$ and
$\mathbf{r}^{\prime}$. To minimize the effect of asymmetry we decided to use a
symmetrized GF

${G}(\mathbf{r},\mathbf{r}^{\prime};\beta)\equiv\frac{1}{2}\left[G(\mathbf{r},\mathbf{r}^{\prime};\beta)+G(\mathbf{r}^{\prime},\mathbf{r};\beta)\right]$
(18)

as our fundamental GF where $G(\mathbf{r},\mathbf{r}^{\prime};\beta)$ in the
RHS of the above equation is defined as

$G(\mathbf{r},\mathbf{r}^{\prime};\beta)={G^{H}}(\mathbf{r},\mathbf{r}^{\prime};\beta)+\int\text{d}\mathbf{r}^{\prime\prime}\bra{\mathbf{r}}e^{-\beta\hat{\mathcal{H}}_{H}}\ket{\mathbf{r}^{\prime\prime}}\bra{\mathbf{r}^{\prime\prime}}\sum_{n=0}^{\infty}\frac{\beta^{n}}{n!}\hat{\mathcal{H}}_{0}^{n}\left(e^{-\beta(\hat{\mathcal{H}}_{0}+\hat{W})}-e^{-\beta\hat{\mathcal{H}}_{0}}\right)\ket{\mathbf{r}^{\prime}}.$
(19)

Eqs.(18) and (19), taken together, are the governing equations to
systematically calculate the KED. The infinite sum over $n$ in Eq.(19)
represents a convergent series of terms thereby facilitating a systematic
scheme of approximations. This feature of our formalism allows for a rigorous
development scheme without re-coursing to any parameter tuning. Each
successive approximation will provide greater accuracy in the similar spirit
of “Jacob’s ladder” in XC functionals. In our scheme, however, no empirical
parameter would be required.

###  II.2 Connection to traditional Green’s function

In GFF we have used Eq.(8) as the basis of our methodology. This form of
Green’s function can be shown to be related to more traditional version of
Green’s functionEconomou (2006) defined as

$\tilde{G}(\mathbf{r},\mathbf{r}^{\prime};E)=\bra{\mathbf{r}}\frac{1}{\hat{\mathcal{H}}-E}\ket{\mathbf{r}^{\prime}}.$
(20)

Numerous methods have been developed to compute $\tilde{G}$ in the problems of
quantum field theory, statistical mechanics, and related fields. As a result a
mathematical connection between $G$ and $\tilde{G}$ will allow for direct
comparison as well as the prospect of employing hitherto existing methods in
our formalism. For example, $\tilde{G}$ can be expanded in terms of the Dyson
equation with the well-understood definition of self-energy $\Sigma$. Highly
effective techniques based on Feynman diagrams have been developed to compute
$\Sigma$. The connection can be achieved by employing an auxiliary variable
$\mu$ _via_

$G(\mathbf{r},\mathbf{r}^{\prime};\beta)=\bra{\mathbf{r}}e^{-\beta\hat{\mathcal{H}}}\ket{\mathbf{r}^{\prime}}=\lim_{\mu\to
1^{+}}\bra{\mathbf{r}}e^{-\mu\beta\hat{\mathcal{H}}}\theta(\mu)\ket{\mathbf{r}^{\prime}}$
(21)

where $\theta(\mu)$ is Heaviside theta function. A subsequent introduction of
Fourier transform $\hat{\mathcal{F}}_{\mu}^{\omega}$ over $\mu$ produces (see
Appendix A for details) the desired relation:

$G(\mathbf{r},\mathbf{r}^{\prime};\beta)=\lim_{\mu\to
1^{+}}\frac{1}{\beta}\hat{\mathcal{F}}_{\mu}^{\omega}(\tilde{G}(\mathbf{r},\mathbf{r}^{\prime};\tau)).$
(22)

The Fourier operator $\hat{\mathcal{F}}_{\mu}^{\omega}$ converts $\mu$ space
functions to $\omega$ space functions and ${\tau=-i\frac{\omega}{\beta}}$.
Since the $\lim_{\mu\to 1^{+}}{}$ and Fourier transform are both linear
operators once can expand $\tilde{G}(\mathbf{r},\mathbf{r}^{\prime};\tau)$
perturbatively before employing $\lim_{\mu\to 1^{+}}$ and
$\hat{\mathcal{F}}_{\mu}^{\omega}$, keeping their relative ordering intact.
The operator $\hat{\tilde{G}}$ can be expanded through Dyson equation

$\hat{\tilde{G}}=\hat{\mathcal{G}}_{0}+\hat{\mathcal{G}}_{0}\hat{W}\hat{\mathcal{G}}_{0}+\hat{\mathcal{G}}_{0}\hat{W}\hat{\mathcal{G}}_{0}\hat{W}\hat{\mathcal{G}}_{0}+\ldots.$
(23)

Here $\hat{\mathcal{G}}_{0}$ is the GF operator for the unperturbed
Hamiltonian which yields
$\mathcal{G}_{0}(\mathbf{r},\mathbf{r}^{\prime};\tau)$ in position
representation and $\hat{W}$ is the mean field interaction potential. The
unperturbed Hamiltonian is $\hat{H}_{H}$ in our case prompting
$\mathcal{G}_{0}(\mathbf{r},\mathbf{r}^{\prime};\tau)\equiv\mathcal{G}^{H}(\mathbf{r},\mathbf{r}^{\prime};\tau)$.
The integral form of the Dyson equation can be written as

$\tilde{G}(\mathbf{r},\mathbf{r}^{\prime};\tau)=\mathcal{G}^{H}(\mathbf{r},\mathbf{r}^{\prime};\tau)+\int\int
d\mathbf{r}^{\prime\prime}d\mathbf{r}^{\prime\prime\prime}\mathcal{G}^{H}(\mathbf{r},\mathbf{r}^{\prime\prime};\tau)W(\mathbf{r}^{\prime\prime},\mathbf{r}^{\prime\prime\prime})\tilde{G}(\mathbf{r}^{\prime\prime\prime},\mathbf{r}^{\prime},\tau).$
(24)

In this work, $\hat{W}$ is considered as a mean field interaction potential
prompting
$W(\mathbf{r}^{\prime\prime},\mathbf{r}^{\prime\prime\prime})\equiv\delta(\mathbf{r}^{\prime\prime}-\mathbf{r}^{\prime\prime\prime})W(\mathbf{r}^{\prime\prime})$.
As a result, Eq.(24) becomes

$\tilde{G}(\mathbf{r},\mathbf{r}^{\prime};\tau)=\mathcal{G}^{H}(\mathbf{r},\mathbf{r}^{\prime};\tau)+\int
d\mathbf{r}^{\prime\prime}\mathcal{G}^{H}(\mathbf{r},\mathbf{r}^{\prime\prime};\tau)W(\mathbf{r}^{\prime\prime})\tilde{G}(\mathbf{r}^{\prime\prime},\mathbf{r}^{\prime},\tau).$
(25)

The expansion for $\tilde{G}$ also has $\mathcal{G}^{H}$ as the leading order
term similar to our GFF while the second term represents a parallel to
$G^{int}$.

##  III Results and Discussion

###  III.1 Weak interaction limit for one-dimensional simple harmonic
oscillator

To exemplify the efficacy of our GFF we apply it to an one-dimensional simple
harmonic oscillator. Also, we have refrained from using the atomic units in
the following sections to clarify the relative magnitudes of different terms.
The ground state of this system with natural frequency $\omega$ and mass $m$
is described by the probability density

$\rho^{SHO}(x)=\left(\frac{m\omega}{\pi\hslash}\right)^{1/2}e^{-\frac{m\omega
x^{2}}{\hslash}}$ (26)

where $\omega$ represents the inverse width of the probability distribution.
As a result, larger $\omega$ leads to a faster decaying probability density
from the center _i.e._ $x=0$. The corresponding GF is

$G(x,x^{\prime};\beta)={G}^{SHO}(x,x^{\prime};\beta)+\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)\bra{x^{\prime\prime}}\sum_{n=0}^{\infty}\frac{\beta^{n}}{n!}\hat{\mathcal{H}}_{0}^{n}\left(e^{-\beta(\hat{\mathcal{H}}_{0}+\hat{W})}-e^{-\beta\hat{\mathcal{H}}_{0}}\right)\ket{x^{\prime}}.$
(27)

In this equation $\hat{W}$ is the mean-field interaction potential. Note that
the nuclear charge ($Z$) in an atom signifies the extent of force that pulls
electrons towards the nucleus in the center. For a harmonic oscillator, a
similar role is played by the force constant $k=m\omega^{2}$. Therefore
Eq.(19) and Eq.(29) are qualitatively equivalent except that for the latter
only bound states are available. Furthermore, we have considered the weak-
field limit since it resembles the low-density limit of the Gell-Mann-Bruckner
modelGell-Mann and Brueckner (1957) where particles move under a constant
positive background leading to constant charge density. As a result, the mean-
field interaction potential also varies negligibly over space. To accomplish
this, we have expanded $\hat{W}$ in power series

$W(x)=W_{0}+x\partial_{x}W(x)|_{0}+\ldots$ (28)

and considered only the first constant term $W_{0}(x)$. This approximation is
equivalent to $\hat{\mathcal{H}}_{0}\approx\hat{\mathcal{H}}_{H}$. By virtue
of this approximation, Eq.(27) becomes

$G(x,x^{\prime};\beta)={G}^{SHO}(x,x^{\prime};\beta)+\int
dx^{\prime\prime}(e^{-\beta
W_{0}}-1)\left({G}^{SHO}(x,x^{\prime\prime};\beta)e^{-\beta\frac{\partial}{\partial\beta}}{G}^{0}(x^{\prime\prime},x^{\prime};\beta)\right).$
(29)

The ${G}^{SHO}(x,x^{\prime\prime};\beta)$ can be computed from path integral
kernel for SHO Feynman  _et al._ (2010) by setting $T=\beta\hslash/i$ to
obtain

${G}^{SHO}(x,x^{\prime\prime};\beta)=A\exp{-B\left[(x^{2}+x^{\prime\prime
2})C-2xx^{\prime\prime}\right]}.$ (30)

Here $A=\sqrt{\frac{m\omega}{2\pi\hslash\sinh({\omega\beta\hslash})}}$,
$B=\frac{m\omega}{2\hslash\sinh({\omega\beta\hslash})}$ and
$C=\cosh({\omega\beta\hslash}).$ The ${G}^{0}(x,x^{\prime\prime};\beta)$ is
similarly obtained as

${G}^{0}(x^{\prime\prime},x^{\prime};\beta)=A_{0}\exp{-B_{0}(x^{\prime\prime}-x^{\prime})^{2}}.$
(31)

Here, $A_{0}=\sqrt{\frac{m}{2\pi\beta\hslash^{2}}}$,
$B_{0}=\frac{m}{2\beta\hslash^{2}}$.

Following the steps elaborated in Appendix B we finally obtain the total GF as

$\bar{G}(x,x^{\prime};\beta)=G^{SHO}(x,x^{\prime};\beta)+\frac{1}{2}\left(\left(\delta-\zeta\right)e^{-\alpha(ax^{2}+ex^{\prime
2}-bxx^{\prime})}+\left(\delta-\zeta^{\prime}\right)e^{-\alpha(ax^{\prime
2}+ex^{2}-bxx^{\prime})}\right)$ (32)

Here

$\displaystyle\delta$ $\displaystyle=\frac{3\left(e^{-\beta
W_{0}}-1\right)}{2}AA_{0}\sqrt{\frac{\pi}{\phi}},$ (33a)
$\displaystyle\zeta(x,x^{\prime})$ $\displaystyle=AA_{0}B_{0}\left(e^{-\beta
W_{0}}-1\right)\sqrt{\frac{\pi}{\phi}}\left[\frac{1}{2\phi}+(\theta(x,x^{\prime}))^{2}\right],$
(33b) $\displaystyle\zeta^{\prime}(x^{\prime},x)$
$\displaystyle=AA_{0}B_{0}\left(e^{-\beta
W_{0}}-1\right)\sqrt{\frac{\pi}{\phi}}\left[\frac{1}{2\phi}+(\theta^{\prime}(x^{\prime},x))^{2}\right],$
(33c) $\displaystyle\phi$ $\displaystyle=(BC+B_{0}),$ (33d)
$\displaystyle\theta(x,x^{\prime})$
$\displaystyle=\frac{B(x-Cx^{\prime})}{\phi},$ (33e)
$\displaystyle\theta^{\prime}(x^{\prime},x)$
$\displaystyle=\frac{B(x^{\prime}-Cx)}{\phi},$ (33f) $\displaystyle a$
$\displaystyle=B^{2}C^{2}-B^{2}+BCB_{0},$ (33g) $\displaystyle e$
$\displaystyle=BCB_{0},$ (33h) $\displaystyle b$ $\displaystyle=2BB_{0},$
(33i) $\displaystyle\alpha$ $\displaystyle=1/\phi.$ (33j)

Following our formalism, we finally obtain the kinetic energy

$T[\rho]=\lim_{\Gamma\to\infty}\frac{1}{2\pi
i}\int_{\gamma-i\Gamma}^{\gamma+i\Gamma}\frac{\text{d}\beta}{\beta}e^{\beta\epsilon_{F}}t(\beta)$
(34)

where $t(\beta)$ is written as

$t(\beta)=\left[\frac{m}{8\omega\hslash^{3}\beta^{2}}+\frac{m\omega}{8\hslash}+\frac{e^{-\beta
W_{0}}-1}{2}\left(\frac{9m\omega}{32\sqrt{2}\hslash}+\frac{3m}{4\sqrt{2}\omega\hslash^{3}\beta^{2}}+\frac{m}{2\sqrt{2}\omega\hslash^{3}\beta^{2}}+\frac{43m\omega}{128\sqrt{2}\hslash}\right)\right]$
(35)

leading to the final KE

$T[\rho]=\frac{m\omega}{8\hslash}+\frac{m\epsilon_{F}^{2}}{4\omega\hslash^{3}}-\frac{5m\epsilon_{F}^{2}}{16\sqrt{2}\omega\hslash^{3}}+\frac{5m(\epsilon_{F}-W_{0})^{2}}{16\sqrt{2}\omega\hslash^{3}}$
(36)

Here $\epsilon_{F}$ is the Fermi energy. We can see clearly from Eq.(36) that
the kinetic energy increases quadratically with
$\epsilon_{F}=\frac{p^{2}_{F}}{2m}$ where $p_{F}$ is the Fermi momentum.
$p_{F}$ signifies the highest possible momentum of the particles in the
system. As a result, the involvement of higher momentum states leads to higher
kinetic energy for the system. On the contrary, the repulsive interaction term
$W_{0}$ pushes particles away from each other as well as from the center of
the SHO. Classically, the momentum of the SHO is maximum at the center. This
effect leads to a lower KE with the increase in $W_{0}$. The variation of the
KE with $\omega$ is more complex. It is evident from Fig.1 that there exists a
frequency $\omega_{min}$ at which the kinetic energy is minimum for a given
$\epsilon_{F}$ and $W_{0}$. It is important to note here that
$W_{0}<<\epsilon_{F}$ should be considered as the weak field limit (see
Sec.III.6 for a more mathematically motivated rationale). Physically, this
definition of weak-field limit can be justified as follows. The Fermi energy
$\epsilon_{F}$ represents the highest energy states for the non-interacting
Hamiltonian which is compensated by central potential provided via a positive
force constant. In realistic atoms, this central potential on negatively
charged electrons is exerted by the positively charged nucleus via Coulomb
interaction. If the inter-particle repulsive interaction $W_{0}$ supersedes
the Fermi energy, the system will become unstable and unbound. As a result,
inter-particle interaction cannot be computed as a perturbative effect. We can
also obtain the expression of $\omega_{min}$ analytically (see Appendix E) as

$\omega_{min}\approx\left(\frac{\epsilon_{F}}{\hslash}\right)\left(2-\frac{5W_{0}}{\sqrt{2}\epsilon_{F}}\right)^{1/2}$
(37)

In other words, the kinetic energy, and by virtue of the virial theorem the
total energy, gets minimized in the presence of interelectronic interaction
for a critical frequency $\omega_{min}\propto\sqrt{k_{min}}$ . The existence
of $\omega_{min}$ is reminiscent of Slater’s shielding constants where the
interelectronic interaction reduces the effective nuclear charges experienced
by an electron in a single-particle state. Since the force constant $k$, and
therefore $\omega$, takes the role of the nuclear charge in our present work
the $\omega_{min}$ hints at the similar phenomenology. Also, as evident from
Eq.(37) this $\omega_{min}$ increases with the $\epsilon_{F}$. This effect is
a parallel to the fact that the shielding due to inner electrons is more
significant than those in higher energy states since these higher energy
electrons are localized farther away from the center of attraction leading to
lesser shielding.

Figure 1: Variation of kinetic energy with $\omega$ for different Fermi
energies $\epsilon_{F}$, keeping $m=1.0$ and $W_{0}=1$. Note that the
$\omega_{\text{min}}$ shifts to higher values with an increase in
$\epsilon_{F}$. The ordinate is plotted on a logarithmic scale.

###  III.2 Application of traditional approach for 1D SHO

A direct term-by-term comparison between our approach and the traditional
Green’s function defined in Eq.(20) is warranted here. Considering the same
approximation of 1D SHO for traditional GF, Eq. 39 can be written as,

$\tilde{G}(x,x^{\prime};\tau)=\mathcal{G}^{H}(x,x^{\prime};\tau)+\int
dx^{\prime\prime}\mathcal{G}^{H}(x,x^{\prime\prime};\tau)W(x^{\prime\prime})\tilde{G}(x^{\prime\prime},x^{\prime},\tau).$
(38)

The idea is to compare our GFF approach with the traditional GF approach using
the same approximation, whereas explained above we replace
$G^{H}(x,x^{\prime};\beta)$ by $G^{SHO}(x,x^{\prime};\beta)$ and
$W(x^{\prime\prime})$ is the mean-field potential as explained in Eq.28, here
also considered the constant interaction potential only.

$\tilde{G}(x,x^{\prime};\tau)=\mathcal{G}^{SHO}(x,x^{\prime};\tau)+\int
dx^{\prime\prime}\mathcal{G}^{SHO}(x,x^{\prime\prime};\tau)W_{0}\tilde{G}(x^{\prime\prime},x^{\prime},\tau).$
(39)

The form of $\mathcal{G}^{SHO}$ with energy $E$ Chua  _et al._ (2018) is

$\mathcal{G}^{SHO}(x,x^{\prime};E)=\frac{m\Gamma(\frac{1}{2}-\epsilon)}{2^{\epsilon-1/2}\sqrt{\pi}\alpha\hslash^{2}}e^{-\frac{\alpha^{2}(x^{2}+x^{\prime
2})}{2}}H_{\epsilon-1/2}(-\alpha\min(x,x^{\prime}))H_{\epsilon-1/2}(\alpha\max(x,x^{\prime})).$
(40)

Here $\alpha=\sqrt{\frac{m\omega}{\hslash}}$,
$\epsilon=\frac{mE}{\alpha^{2}\hslash^{2}}=\frac{E}{\hslash\omega}$ and
$H_{\epsilon-1/2}(x)$ is Hermite polynomial of degree $\epsilon-1/2$. Note
that to compute $\mathcal{G}^{SHO}$, $E$ should be replaced by $\tau$. We
decided to compute the kinetic energy from this expression for the future due
to its complex form. Nevertheless, it is clear from Eq.(40) that the leading
order term of the interaction part of Green’s function varies linearly with
$W_{0}$. This is in contrast with our approach where the leading order
interaction term (See Eq.(33)) is proportional to $\left(e^{-\beta
W_{0}}-1\right)$ thereby incorporating $W_{0}$ contribution up to infinite
order. Also, the $W_{0}$ appears in the exponent ensuring the convergence of
the resulting series. Clearly, each term of our perturbation series appears to
be a partial sum of infinite terms of the Dyson series presented in Eq.(23).
The exact nature of this partial sum, however, requires a sophisticated
Feynman diagram approach which will be taken up in the future.

###  III.3 Spatial nature of Green’s functions

The fundamental entity in our formalism is the Green’s function. Eq.(29) is
the governing equation for GF of a one-dimensional simple harmonic oscillator
(SHO) in the weak interaction limit. In a broad sense, any Green’s function
captures the non-local response of the system, thereby encapsulating all
mechanical information about the system. One can interpret the GF as the
magnitude of the “effect” experienced by an observer at $x^{\prime}$ when a
“cause” has been introduced at $x$. Note that these “cause-effect” points are
arbitrary and therefore the GF is symmetric with respect to $x$ and
$x^{\prime}$. Therefore it is imperative to analyze the spatial nature of the
GF to understand and propel our methodology further. We have plotted Eq.(32)
in Fig.2 for three representative $\omega=0.5,1\text{ and }2$ keeping
$\beta=10$ and $W_{0}=1$ fixed. The contour plots in Fig.2 for $G^{SHO}$,
$G^{int}$, and the total $G$ are generated using different parameter values.
As expected from Eq.(32), all contour plots exhibit $x-x^{\prime}$ symmetry.
Graphically, all contours are symmetrically distributed around the diagonal
line for which $x=x^{\prime}$. The contour of $G^{SHO}$ is completely circular
owing to the fact that the $G^{SHO}(x,x^{\prime},;\beta)$ is a function of
$(x-x^{\prime})$. It also exhibits a maximum response at the center
($x=x^{\prime}=0$) due to the accumulation of probability density
$\rho^{SHO}(x)$ at the center as described by Eq.(26). The response decreases
as one moves away from the center since $\rho^{SHO}(x)$ also tapers off from
the center monotonically. Similarly, $G^{int}$ shows the maximum response at
the center albeit with a negative sign. As a result, the total Green’s
function represents a “tug-of-war” between the unperturbed SHO response and
that of the interaction term. $G^{int}(x,x^{\prime};\beta)$ also tapers off
away from the center due to a decrease in probability density farther away. In
contrast to $G^{SHO}$, $G^{int}$ shows a very prominent “polar” structure. One
can see that with increase in $\omega$ the $G^{int}$ becomes more localized
mirroring $\rho^{SHO}(x)$. Clearly, the response is longer in range along
$x=0$ and $x^{\prime}=0$ lines. This feature becomes more prominent with the
increase in $\omega$. While the exact reason for this feature is not clear at
present, we believe that it is related to the density difference between
“cause” and “effect” points. Since $x=0$ or $x^{\prime}=0$, line represents
the region maximum probability of the particle, $G^{int}$ along these lines
indicate a large difference of density leading to stronger interaction
indicating a strong correlation between the density gradient
$\frac{\partial\rho^{SHO}(x)}{\partial x}$ and interparticle interaction
$W_{0}$. This feature is clearly reminiscent of inhomogeneity in atomic
electron density exemplified by atomic shell structures. Clearly, this
interaction term appears to hold the key to atomic shell structures. However,
more work is necessary to prove this conclusion explicitly. The total $G$ is a
combined effect of $G^{SHO}$ and $G^{int}$ and is more difficult to understand
physically. It is clear, however, that while the $G^{int}$ becomes less
“polar” with the decrease in $\omega$, total GF $G(x,x^{\prime};\beta)$
becomes more “polar” indicating a complex interplay between $G^{SHO}$ and
$G^{int}$. Future works are essential to understand these features for
devising novel methodologies applicable to realistic systems.

|  A |  B |  C
---|---|---|---
$\mathbf{\omega=2}$ |  |  |
$\mathbf{\omega=1}$ |  |  |
$\mathbf{\omega=0.5}$ |  |  |
Figure 2: The contour plots of $G^{\text{SHO}}$, $G^{\text{int}}$, and $G$ for
different $\omega$ values and a constant $\beta=10$. The first column (A)
displays $G^{\text{SHO}}$, the second column (B) shows $G^{\text{int}}$, and
the third column (C) illustrates $G$. The rows correspond to different
$\omega$ values: the first row is for $\omega=2$, the second row is for
$\omega=1$, and the third row is for $\omega=0.5$. Each contour plot provides
a visual representation of the corresponding GFs under the specified
parameters.

###  III.4 $N$ and $Z$ dependence of KE

Next, we will discuss the variation of the KE expression obtained in Eq.(36)
with the number of particle $N$ via $\epsilon_{F}$ and $\omega$. We realize
that $\epsilon_{F}$ indicates the upper limit of the occupied energy states
and is therefore indicative of $N$ in a given set of single-particle states.
We have employed two estimates of $\epsilon_{F}$ in this work resulting in two
asymptotic limits. They represent pure bound state limit (BSL) and pure
scattering state limit (SSL). This is done keeping in mind that realistic
systems like hydrogenic atoms sustain both types of states and therefore these
two models act as as asymptotic limits of the real atoms. Also, we considered
each single-particle state to be singly occupied to mimic the effects of the
Pauli exclusion principle. First, we consider the quantum harmonic oscillator
as a representative of BSL to define $\epsilon_{F}$. As a result, we define
bound state limit (BSL) by $\epsilon_{F}=(N+1/2)\hslash\omega$. The kinetic
energy functional for BSL is therefore

$T_{BSL}\approx A(N+1/2)^{2}-B(N+1/2)+C+D$ (41)

where $m$ and $\omega$ dependent constants A ,B, C and D are

$\displaystyle A$ $\displaystyle=$ $\displaystyle\frac{m\omega}{4\hslash}$
(42) $\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{5mW_{0}}{8\sqrt{2}\hslash^{2}}$ (43) $\displaystyle C$
$\displaystyle=$ $\displaystyle\frac{m\omega}{8\hslash}$ (44) $\displaystyle
D$ $\displaystyle=$
$\displaystyle\frac{5mW_{0}^{2}}{16\sqrt{2}\omega\hslash^{3}}.$ (45)

The SSL is represented by free particles with a continuous energy spectrum. In
this limit

$\epsilon_{F}=\frac{\hslash^{2}}{2m}\left(\frac{3\pi^{2}}{V}\right)^{2/3}N^{2/3}$
(46)

where the volume $V$ is computed here via classical turning points $x_{tp}$ of
the SHO i.e $V=(2x_{tp})^{3}$. We employed the classical turning points are
related to the energy of the highest occupied state
$E_{N}=\left(N+\frac{1}{2}\right)\hslash\omega$ and given by

$x_{tp}=\sqrt{\frac{2E_{N}}{m\omega^{2}}}$ (47)

As a result, we obtain

$V=8\left(\frac{2(N+1/2)\hslash}{m\omega}\right)^{3/2}$ (48)

leading to the $\epsilon_{F}$

$\epsilon_{F}=\frac{3\pi^{2}\hslash\omega}{16}\frac{N^{2/3}}{(N+1/2)}.$ (49)

Employing these expressions we finally find

$T_{SSL}\approx E\frac{N^{4/3}}{(N+1/2)^{2}}-F\frac{N^{2/3}}{(N+1/2)}+C+D.$
(50)

where

$\displaystyle E$ $\displaystyle=$
$\displaystyle\left(\frac{3\pi^{2}}{16}\right)^{2}\left(\frac{m\omega}{4\hslash}\right)$
(51) $\displaystyle F$ $\displaystyle=$
$\displaystyle\left(\frac{3\pi^{2}}{16}\right)\left(\frac{5mW_{0}}{8\sqrt{2}\hslash^{2}}\right)$
(52)

To compare the relative qualitative accuracy of the four aforementioned KE
functionals, we employed empirical $N$ and $Z$-dependency of kinetic energy
$T[\rho]=T_{vW}[\rho]+\gamma(N,Z)T_{0}[\rho]$ defined for atomic systems. For
neutral atoms $N=Z$ leads to $\gamma(N,N)=1-\frac{1.412}{N^{1/3}}$Parr and
Yang (1994); Acharya  _et al._ (1980). Furthermore, to get the $N-$dependence
of KE we employ $\rho\propto N$ to finally obtain

$T_{\gamma}\approx N^{2}+N^{5/3}-(1.412)N^{4/3}.$ (53)

Figure 3: The variation of kinetic energy as a function of $N$. The blue solid
line, yellow solid line, and black dashed line represent $T_{SSL}$, $T_{BSL}$
and $T_{\gamma}$, respectively. The ordinate is in logscale.

A direct comparison between $T_{BSL}$, $T_{SSL}$ and $T_{\gamma}$ in Fig.3. It
is evident that the $T_{BSL}$ follows the same functional dependence on $N$ as
$T_{\gamma}$. The behavior of $T_{SSL}$ is completely wrong. This is expected
since the realistic atoms are closer in essence to BSL than SSL. The
difference between $T_{\gamma}$ and $T_{BSL}$ also increases for larger $N$.
This effect is also expected since the single-particle states in SHO remain
equispaced in energy while for atoms they become more closely spaced.
Therefore for larger $N$ the qualitative difference between our model and
actual atoms also gets magnified leading to larger error. While a quantitative
agreement between our model and an actual atom is completely unexpected, our
model appears to reproduce the qualitative features of realistic atoms,
including the $N$-dependence of the KE, quite successfully.

###  III.5 Application of $T_{BSL}$ on noble gas atoms

To push our formalism further, we have now applied our results to noble gas
atoms. Note that it is expected that our results will not match the actual
kinetic energy numerically for our expression is derived for a one-dimensional
(1D) system. and incorporates other crude approximations. However, we are
interested in finding out if there is any qualitative trend that emerges from
$(T_{BSL})$ applied to noble gas atoms using Eq.(41). The reference kinetic
energy values were taken from Hartree-Fock table data Clementi and Roetti
(1974). We employed two estimates for $\omega$. First, we have used the
ionization energies of the corresponding atoms in atomic units (a.u.) as an
estimate of $\omega=\omega_{IE}$ Kramida  _et al._ (2023). The second estimate
of $\omega=\omega_{K-L}$ is taken to be the transition frequency between K and
L shells of the respected atoms. Results from these two estimates are
presented in Table.1. Note that except for He where these are essentially the
same due to following Koopmans theorem, for other atoms
$\omega_{IE}<\omega_{K-L}$. Moreover $\omega_{IE}$ and $\omega_{K-L}$
represent the smallest and largest energy gaps for these atoms. Therefore it
is expected that actual values of the KE will lie between these two extreme
limits. This expectation is indeed found to be true. It can be seen from Table
1 that $T_{BSL}$ and $T_{BSL}^{K-L}$ represent lower and higher extremes of
true KE. In future, a search for an effective frequency may be pursued which
will reproduce the KE more faithfully. Such quantity will turn out to be very
useful for Unsold-type approximations in the computations of effective
response functions Wiles and Manby (2018).

Atom |  $\mathbf{\omega_{IE}}$ (a.u.) |  $\mathbf{T_{BSL}}$ (a.u.) |  $\mathbf{\omega_{K-L}}$ (a.u.) |  $\mathbf{T_{BSL}^{K-L}}$ (a.u.) |  K.E (a.u.)
---|---|---|---|---|---
He | $9.1\times 10^{-1}$ | $6.7\times 10^{-1}$ | $9.2\times 10^{-1}$ | $6.8\times 10^{-1}$ | $2.86\times 10^{0}$
Ne | $7.9\times 10^{-1}$ | $1.76\times 10$ | $3.08\times 10^{0}$ | $8.49\times 10^{2}$ | $1.28\times 10^{2}$
Ar | $5.8\times 10^{-1}$ | $4.18\times 10$ | $1.06\times 10^{2}$ | $9.09\times 10^{3}$ | $5.26\times 10^{2}$
Kr | $5.1\times 10^{-1}$ | $1.55\times 10^{2}$ | $4.50\times 10^{2}$ | $1.5\times 10^{5}$ | $2.75\times 10^{3}$
Table 1: Calculated kinetic energy values, $T_{BSL}$ and $T_{BSL}^{K-L}$, for
noble gas atoms using $\omega_{IE}$ and $\omega_{K-L}$, respectively. The
reference kinetic energy (K.E.) values are taken from Hartree-Fock table data
Clementi and Roetti (1974). The parameters $\omega_{IE}$Kramida  _et al._
(2023) and $\omega_{K-L}$ are used to compute the $T_{BSL}$ and
$T_{BSL}^{K-L}$ values, respectively.

###  III.6 A natural limit of weak-field interaction

We have explored the weak field limit extensively in our work and defined it
as

$W_{0}^{max}<\epsilon_{F}$ (54)

in Section III.1 from a physical consideration where $W_{0}^{max}$ is the
maximum allowed value for interaction potential. In addition, our formalism
supports this definition in a more mathematically rigorous manner. The
condition to obtain a well-defined Laplace transform in Eq.(34) is

$\epsilon_{F}-W_{0}\geq 0.$ (55)

Without the condition, the Bromwich integral does not yield a bounded KE
justifying our physical definition in Eq.(54). Note that $\epsilon_{F}$
increases with $N$. Therefore for $N\to\infty$, any system would be considered
as the weak-field limit leading to an exact HEG description. This conclusion
coincides with Lieb’s resultsLieb (1976).

##  IV Conclusion

In summary, we have developed a general and systematic method to compute the
kinetic energy density for an electronic system using Green’s function
formalism. The methodology developed in this work is the first approach, to
our knowledge, to compute KED in a completely systematic manner. Another
significant advantage of our approach is that it is not based on top of an
arbitrary model. Instead, the guiding principle of our perturbative treatment
is the magnitude of the different partitions of the total Hamiltonian leading
to a mathematically consistent methodology for KED. We have also shown that
our methodology is related to traditional Green’s function formalism in a
mathematically rigorous manner. However, each term of our expansion can be
considered as a partial sum of the traditional perturbation series of Green’s
functions. To exemplify the formalism we also presented some flagship
calculations on a one-dimensional Hookian atom. While we did not expect any
quantitative agreement with actual kinetic energy values of atoms these
calculations successfully reproduced several salient qualitative features of
atomic KEDs such as the the variation of KE on $N$. To the best of our
knowledge, this is the first work which shows such dependence without
recoursing to any specific model. From the calculation of KE for weak
interaction limit, we also obtained an upper limit of the magnitude of
interaction potential in terms of the Fermi energy. Such quantitative measure
for “weak” potential is also one of the natural by-products of our formalism.
Furthermore, the mathematically rigorous connection between traditional
Green’s function and our Green’s function allows for an exact recipe for the
systematic development of KE.

Despite being promising, our methodology requires improvements in several
aspects. First, the form of self-energy computed in this work is extremely
complicated and only computed in a mean-field framework. A more rigorous many-
body description needs to be developed in the future. A full Fermionic many-
body expansion of the Green’s function in its present form appears
challenging. Therefore, in the future, more sophisticated approximation
techniques need to be developed in that direction. In the near future, the use
of Hydrogenic Green’s functions also have to be used and research in that
direction is in progress. Moreover, a diagrammatic perturbation technique is
essential to develop an analytically convergent perturbation theory.

###### Acknowledgements.

Priya acknowledges a Prime Minister’s Research Fellowship and MS acknowledges
IIT Kanpur for initiation grant for partial support. We also thank Jaseela for
support.

##  Appendix A Connection to $G(\mathbf{r},\mathbf{r}^{\prime};\beta)$ and
$\tilde{G}(\mathbf{r},\mathbf{r}^{\prime};\tau)$

We define the Fourier transform operator $\hat{\mathcal{F}_{\mu}^{\omega}}$ as

$F(\omega)=\hat{\mathcal{F}_{\mu}^{\omega}}(f(\mu))\equiv\frac{1}{\sqrt{2\pi}}\int
e^{i\mu\omega}f(\mu)\text{d}\mu.$ (56)

Using Eq.(56), we can our Green’s function in terms of traditional GF as

$\displaystyle G(\mathbf{r},\mathbf{r}^{\prime};\beta)$ $\displaystyle=$
$\displaystyle\lim_{\mu\to
1^{+}}\bra{\mathbf{r}}\hat{\mathcal{F}}_{\mu}^{\omega}(\frac{1}{\beta\hat{\mathcal{H}}+i\omega})\ket{\mathbf{r}^{\prime}}$
(57) $\displaystyle=$ $\displaystyle\lim_{\mu\to
1^{+}}\frac{1}{\beta}\bra{\mathbf{r}}\hat{\mathcal{F}}_{\mu}^{\omega}(\frac{1}{\hat{\mathcal{H}}+i\frac{\omega}{\beta}})\ket{\mathbf{r}^{\prime}}$
$\displaystyle=$ $\displaystyle\lim_{\mu\to
1^{+}}\frac{1}{\beta}\hat{\mathcal{F}}_{\mu}^{\omega}(\tilde{G}(\mathbf{r},\mathbf{r}^{\prime};\tau)).$

where $\tau=\frac{-i\omega}{\beta}$.

##  Appendix B Calculation for 1D harmonic oscillator using GFF

The detailed calculation for kinetic energy is described here.

$G(x,x^{\prime};\beta)={G}^{SHO}(x,x^{\prime};\beta)+\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)\bra{x^{\prime\prime}}\sum_{n=0}^{\infty}\frac{\beta^{n}}{n!}\hat{\mathcal{H}}_{0}^{n}\left(e^{-\beta(\hat{\mathcal{H}}_{0}+\hat{W})}-e^{-\beta\hat{\mathcal{H}}_{0}}\right)\ket{x^{\prime}}.$
(58)

Basic definition used in our work

${G}^{SHO}(x,x^{\prime\prime};\beta)=A\exp{-B\left[(x^{2}+x^{\prime\prime
2})C-2xx^{\prime\prime}\right]}.$ (59)

Here

$\displaystyle A$ $\displaystyle=$
$\displaystyle\sqrt{\frac{m\omega}{2\pi\hslash\sinh({\omega\beta\hslash})}}$
$\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{m\omega}{2\hslash\sinh({\omega\beta\hslash})}$
$\displaystyle C$ $\displaystyle=$ $\displaystyle\cosh({\omega\beta\hslash})$

The ${G}^{0}(x,x^{\prime\prime};\beta)$ is similarly obtained as

${G}^{0}(x^{\prime\prime},x^{\prime};\beta)=A_{0}\exp{-B_{0}(x^{\prime\prime}-x^{\prime})^{2}}.$
(60) $\displaystyle A_{0}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{m}{2\pi\beta\hslash^{2}}}$ $\displaystyle B_{0}$
$\displaystyle=$ $\displaystyle\frac{m}{2\beta\hslash^{2}}.$

as explained above in section III, we considered weak interaction limit

$W(x)=W_{0}+x\partial_{x}W(x)|_{0}+\ldots$ (61)

We have considered the constant interaction potential, as a result eq. (27)
becomes

$G(x,x^{\prime};\beta)={G}^{SHO}(x,x^{\prime};\beta)+\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)\bra{x^{\prime\prime}}\sum_{n=0}^{\infty}\frac{\beta^{n}}{n!}\hat{\mathcal{H}}_{0}^{n}\left(e^{-\beta\hat{\mathcal{H}}_{0}}e^{-\beta\hat{W}_{0}}-e^{-\beta\hat{\mathcal{H}}_{0}}\right)\ket{x^{\prime}}.$
(62)

rewriting the above equation

$G(x,x^{\prime};\beta)={G}^{SHO}(x,x^{\prime};\beta)+\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)(e^{-\beta
W_{0}}-1)\bra{x^{\prime\prime}}\sum_{n=0}^{\infty}\frac{\beta^{n}}{n!}\underbrace{\hat{\mathcal{H}}_{0}^{n}e^{-\beta\hat{\mathcal{H}}_{0}}}_{\frac{\partial^{n}}{\partial(-\beta)^{n}}e^{-\beta\hat{\mathcal{H}}_{0}}}\ket{x^{\prime}}.$
(63)

$G(x,x^{\prime};\beta)={G}^{SHO}(x,x^{\prime};\beta)+(e^{-\beta W_{0}}-1)\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)e^{-\beta\frac{\partial}{\partial\beta}}G^{0}(x^{\prime\prime},x^{\prime};\beta)$
(64)

solving for the 2nd part
$e^{-\beta\frac{\partial}{\partial\beta}}G^{0}(x^{\prime\prime},x^{\prime};\beta)\approx\left(1-\beta\frac{\partial}{\partial\beta}\right)G^{0}(x^{\prime\prime},x^{\prime};\beta)=\frac{3}{2}G^{0}(x^{\prime\prime},x^{\prime};\beta)-B_{0}(x^{\prime\prime}-x^{\prime})^{2}G^{0}(x^{\prime\prime},x^{\prime};\beta)$

Therefore we get

$\begin{split}G(x,x^{\prime};\beta)={G}^{SHO}(x,x^{\prime};\beta)+(e^{-\beta
W_{0}}-1)\left(\frac{3}{2}\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)G^{0}(x^{\prime\prime},x^{\prime};\beta)\right)\\\
-(e^{-\beta W_{0}}-1)\left(B_{0}\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)(x^{\prime\prime}-x^{\prime})^{2}G^{0}(x^{\prime\prime},x^{\prime};\beta)\right)\end{split}$
(65)

##  Appendix C The kinetic energy contribution from $G^{SHO}$

The kinetic energy contribution from $G^{SHO}$ in Eq.(65) can be calculated as

$\displaystyle T[\rho]$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int{\frac{\partial^{2}}{\partial
x^{2}}(\rho_{1}(x,x^{\prime}))}|_{x=x^{\prime}}dx$ (66)

As explained in section II $\rho_{1}(x,x^{\prime})$ is related to GF by
Laplace transform

$\displaystyle T[\rho]$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int{\frac{\partial^{2}}{\partial
x^{2}}}\left(\lim_{T\to\infty}\frac{1}{2\pi i}\int_{\gamma-
iT}^{\gamma+iT}\frac{\text{d}\beta}{\beta}e^{\beta\epsilon_{F}}G(x,x^{\prime};\beta)\right)|_{x=x^{\prime}}dx$
(67) $\displaystyle=$ $\displaystyle-\frac{1}{2}\lim_{T\to\infty}\frac{1}{2\pi
i}\int_{\gamma-
iT}^{\gamma+iT}\frac{\text{d}\beta}{\beta}e^{\beta\epsilon_{F}}\int{\frac{\partial^{2}}{\partial
x^{2}}}G(x,x^{\prime};\beta)|_{x=x^{\prime}}dx$ (68)

Here we are doing for $G^{SHO}$, therefore we need to evaluate
${\frac{\partial^{2}}{\partial x^{2}}}G^{SHO}(x,x^{\prime};\beta)$

$\frac{\partial^{2}}{\partial
x^{2}}G^{SHO}(x,x^{\prime};\beta)|_{x=x^{\prime}}=4AB^{2}x^{2}(C-1)^{2}e^{-2Bx^{2}(C-1)}-2ABCe^{-2Bx^{2}(C-1)}$
(69)

$\frac{\partial^{2}}{\partial x^{\prime
2}}G^{SHO}(x,x^{\prime};\beta)=\frac{\partial^{2}}{\partial
x^{2}}G^{SHO}(x,x^{\prime};\beta)$ (70)

Now doing the integration w.r.t x

$\displaystyle\int\frac{\partial^{2}}{\partial
x^{2}}G^{SHO}(x,x^{\prime};\beta)|_{x=x^{\prime}}dx$ $\displaystyle=$
$\displaystyle\int
dx(4AB^{2}x^{2}(C-1)^{2}e^{-2Bx^{2}(C-1)}-2ABCe^{-2Bx^{2}(C-1)})$ (71)
$\displaystyle=$ $\displaystyle
4AB^{2}(C-1)^{2}\frac{1}{2}\sqrt{\frac{\pi}{(2B(C-1))^{3})}}-2ABC\sqrt{\frac{{\pi}}{{2B(C-1)}}}$
(72) $\displaystyle=$ $\displaystyle
A\sqrt{B\pi}\left(\sqrt{\frac{C-1}{2}}-\sqrt{\frac{2}{C-1}}C\right)$ (73)
$\displaystyle=$ $\displaystyle-A\sqrt{B\pi}\frac{C+1}{\sqrt{2(C-1)}}$ (74)

Now using $A=\sqrt{B/\pi}$ , $C=1+\frac{(\omega\beta\hslash)^{2}}{2}$ and
using approximation $sinh(\omega\beta\hslash)=\omega\beta\hslash$, therefore
$B=\frac{m}{2\beta\hslash^{2}}$. Using all the approximations we get

$\displaystyle=$
$\displaystyle-\frac{m}{4\beta\hslash^{2}}\left(\frac{4}{\omega\beta\hslash}+\omega\beta\hslash\right)$
(75) $\displaystyle=$ $\displaystyle\frac{1}{2}\lim_{T\to\infty}\frac{1}{2\pi
i}\int_{\gamma-
iT}^{\gamma+iT}\frac{\text{d}\beta}{\beta}e^{\beta\epsilon_{F}}\frac{m}{4\beta\hslash^{2}}\left(\frac{4}{\omega\beta\hslash}+\omega\beta\hslash\right)$
(76) $\displaystyle=$
$\displaystyle\frac{m}{8\hslash^{2}}\left(\lim_{T\to\infty}\frac{1}{2\pi
i}\int_{\gamma-iT}^{\gamma+iT}\text{d}\beta
e^{\beta\epsilon_{F}}\left(\frac{4}{\omega\beta^{3}\hslash}+\frac{\omega\hslash}{\beta}\right)\right)$
(77) $\displaystyle T[\rho]_{SHO}$ $\displaystyle=$
$\displaystyle\frac{m\epsilon_{F}^{2}}{4\omega\hslash^{3}}+\frac{m\omega}{8\hslash}$
(78)

##  Appendix D The Kinetic energy contribution from $G^{int}$ in Eq.(65)

Another essential required for solving the above equation is the
multiplication of two Gaussian in our case
$G^{SHO}(x,x^{\prime\prime};\beta)G^{0}{x^{\prime\prime},x^{\prime};\beta}$
writing the general form

$\displaystyle e^{-\alpha_{1}(x-\mu_{1})^{2}}e^{-\alpha_{2}(x-\mu_{2})^{2}}$
$\displaystyle=$ $\displaystyle
e^{-(\alpha_{1}\mu_{1}^{2}+\alpha_{2}\mu_{2}^{2})}e^{-\left(\frac{(\alpha_{1}\mu_{1}+\alpha_{2}\mu_{2})^{2}}{\alpha_{1}+\alpha_{2}}\right)}e^{-(\alpha_{1}+\alpha_{2})(x-z)^{2}}$
(79)

Here $z=\frac{(\alpha_{1}\mu_{1}+\alpha_{2}\mu_{2})}{\alpha_{1}+\alpha_{2}}$.
Similarly, we have to solve another general result for 3rd term

$\displaystyle\int
dx^{\prime\prime}(x^{\prime\prime}-\theta_{1})^{2}e^{-\phi(x^{\prime\prime}-\theta_{2})^{2}}$
$\displaystyle=$ $\displaystyle\int
dx^{\prime\prime}((x-\theta_{2})+\theta)^{2}e^{-\phi(x^{\prime\prime}-\theta_{2})^{2}}$
(80) $\displaystyle=$
$\displaystyle\sqrt{\frac{\pi}{\phi}}\left[\frac{1}{2\phi}+\theta^{2}\right]$
(81)

Using eq.(79) the multiplication of
$G^{SHO}(x,x^{\prime\prime};\beta)G^{0}(x^{\prime\prime},x^{\prime}:\beta)$ is

$\displaystyle
G^{SHO}(x,x^{\prime\prime};\beta)G^{0}(x^{\prime\prime},x^{\prime}:\beta)$
$\displaystyle=$ $\displaystyle
e^{-\phi(x^{\prime\prime}-z)^{2}}e^{-\alpha(ax^{2}+ex^{\prime
2}-bxx^{\prime})}$ $\displaystyle\phi$ $\displaystyle=$
$\displaystyle(BC+B_{0})$ $\displaystyle a$ $\displaystyle=$ $\displaystyle
B^{2}C^{2}-B^{2}+BCB_{0}$ $\displaystyle e$ $\displaystyle=$ $\displaystyle
BCB_{0}$ $\displaystyle b$ $\displaystyle=$ $\displaystyle 2BB_{0}$
$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle 1/\phi$ $\displaystyle
z$ $\displaystyle=$ $\displaystyle\frac{Bx+B_{0}x^{\prime}}{BC+B_{0}}$

Solving for the 1st part of $G^{int}$ in Eq.(65)

$\displaystyle\int
dx^{\prime\prime}G^{SHO}(x,x^{\prime\prime};\beta)G^{0}(x^{\prime\prime},x^{\prime}:\beta)$
$\displaystyle=$ $\displaystyle
AA_{0}\sqrt{\frac{\pi}{BC+B_{0}}}e^{-\alpha(ax^{2}+ex^{\prime
2}-bxx^{\prime})}$ Now double differentiating $e^{-\alpha(ax^{2}+ex^{\prime
2}-bxx^{\prime})}$ w.r.t $x$ and $x=x^{\prime}$
$\displaystyle\left(\alpha^{2}(2a-b)^{2}x^{2}-2\alpha
a\right)e^{-\alpha.x^{2}(a+e-b)}$ Now double differentiating
$e^{-\alpha(ax^{2}+ex^{\prime 2}-bxx^{\prime})}$ w.r.t $x^{\prime}$ and
$x=x^{\prime}$ $\displaystyle\left(\alpha^{2}(2e-b)^{2}x^{2}-2\alpha
e\right)e^{-\alpha.x^{2}(a+e-b)}$ Now Taking the average
$\displaystyle\frac{1}{2}\left(\alpha^{2}(2a-b)^{2}x^{2}-2\alpha
a+\alpha^{2}(2e-b)^{2}x^{2}-2\alpha e\right)e^{-\alpha.x^{2}(a+e-b)}$
Integrating w.r.t x
$\displaystyle\frac{\alpha^{2}(2a-b)^{2}\pi^{1/2}}{2(\alpha(a+e-b))^{3/2}}-\frac{2\alpha
a\pi^{1/2}}{(a+e-b)^{1/2}}+\frac{\alpha^{2}(2e-b)^{2}\pi^{1/2}}{2(\alpha(a+e-b))^{3/2}}-\frac{2\alpha
e\pi^{1/2}}{(a+e-b)^{1/2}}$

Including the $AA_{0}\sqrt{\frac{\pi}{BC+B_{0}}}$ and using the approximations
discussed above which makes $A=A_{0},B=B_{0}$ we get

$\frac{3(e^{-\beta
W_{0}}-1)}{4(C+1)(a+e-b)^{1/2}}\left[\frac{(2a-b)^{2}}{2(a+e-b)}+\frac{(2e-b)^{2}}{2(a+e-b)}-2a-2e\right]$
(82)

Using $C=(1+\frac{(\omega\beta\hslash)^{2}}{2})$

$\displaystyle a+e-b$ $\displaystyle=$ $\displaystyle
2B^{2}(\omega\beta\hslash)^{2}$ $\displaystyle 2a-b$ $\displaystyle=$
$\displaystyle 3B^{2}(\omega\beta\hslash)^{2}$ $\displaystyle 2e-b$
$\displaystyle=$ $\displaystyle B^{2}(\omega\beta\hslash)^{2}$ $\displaystyle
2e$ $\displaystyle=$ $\displaystyle 2B^{2}C$ $\displaystyle 2a$
$\displaystyle=$ $\displaystyle
2B^{2}\left(1+\frac{3(\omega\beta\hslash)^{2}}{2}\right)$

By using the above equations, the Eq. 82 becomes

$-\frac{3(e^{-\beta
W_{0}}-1)}{8\sqrt{2}}\frac{m}{\hslash^{2}}\left(\frac{3\omega\hslash}{4}+\frac{2}{\omega\hslash\beta^{2}}\right)$
(83)

Hereafter multiplying by the factor of $1/\beta$ before doing the inverse
Laplace transform

$-\frac{3m}{8\sqrt{2}\hslash^{2}}\left[\frac{3\omega\hslash e^{-\beta
W_{0}}}{4\beta}+\frac{2e^{-\beta
W_{0}}}{(\omega\hslash\beta^{3})}-\frac{3\omega\hslash}{4\beta}-\frac{2}{\omega\hslash\beta^{3}}\right]$
(84)

The kinetic energy expression that we got is

$T[\rho]_{1}=\frac{3m}{16\sqrt{2}\omega\hslash^{3}}\left((\epsilon_{F}-W_{0})^{2}-\epsilon_{F}^{2}\right)$
(85)

####  D.0.1 The kinetic energy contribution from the second part $G^{int}$ of
Eq.(65)

The second part of Eq. (65) is

$=(e^{-\beta W_{0}}-1)\left(B_{0}\int
dx^{\prime\prime}{G}^{SHO}(x,x^{\prime\prime};\beta)(x^{\prime\prime}-x^{\prime})^{2}G^{0}(x^{\prime\prime},x^{\prime};\beta)\right)$
(86)

Using the general formula Eq.(79) and Eq.(80) we get

$\displaystyle\sqrt{\frac{\pi}{\phi}}\left(\frac{1}{2\phi}+\theta^{2}\right)e^{-}\left({\frac{ax^{2}+ex^{\prime
2}-bxx^{\prime}}{\phi}}\right)$
$\displaystyle\theta=\frac{B(x-Cx^{\prime})}{\phi}$ Simplifying and
substituting the value of $\theta$
$\displaystyle\left(L+M(x-Cx^{\prime})\right)e^{-}\left({\frac{ax^{2}+ex^{\prime
2}-bxx^{\prime}}{\phi}}\right)$ here $L$ and $M$ are $\displaystyle
L=\frac{\pi^{1/2}}{2\phi^{3/2}};M=\frac{B^{2}\pi^{1/2}}{\phi^{5/2}}$

Now further we have to double differentiate w.r.t $x$ and $x^{\prime}$

$(L+M(x-Cx^{\prime})^{2})e^{-}\left({\frac{ax^{2}+ex^{\prime
2}-bxx^{\prime}}{\phi}}\right)$ (87)

After differentiating w.r.t x

$[2M-4\alpha M(x-Cx^{\prime})(2ax-bx^{\prime})-2\alpha
a(L+M(x-Cx^{\prime})^{2})+\alpha^{2}(L+M(x-Cx^{\prime})^{2})(2ax-
bx^{\prime})^{2}]e^{-\alpha(ax^{2}+ex^{\prime 2}-bxx^{\prime})}$ (88)

Similarly differentiating w.r.t x’

$[2MC^{2}+4MC\alpha(x-Cx^{\prime})(2ex^{\prime}-bx)-2\alpha
e(L+M(x-Cx^{\prime})^{2})+\alpha^{2}(2ex^{\prime}-bx)^{2}(L+M(x-Cx^{\prime})^{2})]e^{-\alpha(ax^{2}+ex^{\prime
2}-bxx^{\prime})}$ (89)

After using the approximations $sinh(\omega\hslash\beta)=\omega\hslash\beta$,
$cosh(\omega\hslash\beta)=1+\frac{(\omega\hslash\beta)^{2}}{2}$ and
subsituting the values of $L,M,C,\alpha$ etc.. Putting them together in
Eq.(86) we get

$(e^{-\beta
W_{0}}-1)\left(\frac{m}{16\sqrt{2}\omega\hslash^{3}}\left(\frac{8}{\beta^{2}}+\frac{43(\omega\hslash)^{2}}{8}\right)\right)$
(90)

After doing the Laplace transform

$T[\rho]_{2}=\frac{m}{8\sqrt{2}\omega\hslash^{3}}\left((\epsilon_{F}-W_{0})^{2}-\epsilon_{F}^{2}\right)$
(91)

##  Appendix E Calculation of $\omega_{min}$

In order to calculate the $\omega_{min}$ form the kinetic energy (Eq.(36))

$T[\rho]=\frac{m\omega}{8\hslash}+\frac{m\epsilon_{F}^{2}}{4\omega\hslash^{3}}-\frac{5m\epsilon_{F}^{2}}{16\sqrt{2}\omega\hslash^{3}}+\frac{5m(\epsilon_{F}-W_{0})^{2}}{16\sqrt{2}\omega\hslash^{3}}$
(92)

Minimize the $T[\rho]$ w.r.t $\omega$

$\frac{\partial
T[\rho]}{\partial\omega}=\frac{m}{8\hslash}-\frac{m\epsilon_{F}^{2}}{4\omega^{2}\hslash^{3}}+\frac{5m\epsilon_{F}^{2}}{16\sqrt{2}\omega^{2}\hslash^{3}}-\frac{5m(\epsilon_{F}-W_{0})^{2}}{16\sqrt{2}\omega^{2}\hslash^{3}}=0$
(93)

The resulting expression for $\omega_{min}$ is

$\omega_{min}=\frac{\epsilon_{F}}{\hslash}\left(2-\frac{5W_{0}}{\sqrt{2}\epsilon_{F}}\right)^{1/2}$
(94)

## References

  * (1) E. V. Ludeña and V. V. Karasiev, “Kinetic energy functionals: History, challenges and prospects,” in  _Reviews of Modern Quantum Chemistry_, pp. 612–665.
  * Wesolowski and Wang (2013) T. A. Wesolowski and Y. A. Wang,  _Recent Progress in Orbital-free Density Functional Theory_ (WORLD SCIENTIFIC, 2013) https://www.worldscientific.com/doi/pdf/10.1142/8633 .
  * García-Aldea and Alvarellos (2008) D. García-Aldea and J. E. Alvarellos, Phys. Rev. A 77, 022502 (2008).
  * Perdew and Constantin (2007) J. P. Perdew and L. A. Constantin, Phys. Rev. B 75, 155109 (2007).
  * Ghiringhelli and Delle Site (2008) L. M. Ghiringhelli and L. Delle Site, Phys. Rev. B 77, 073104 (2008).
  * Eek and Nordholm (2006) W. Eek and S. Nordholm, Theoretical Chemistry Accounts 115, 266 – 273 (2006), cited by: 10.
  * Karasiev  _et al._ (2006) V. Karasiev, S. Trickey, and F. E. Harris, Journal of Computer-Aided Materials Design 13, 111 – 129 (2006), cited by: 55.
  * Karasiev  _et al._ (2009) V. V. Karasiev, R. S. Jones, S. B. Trickey, and F. E. Harris, Phys. Rev. B 80, 245120 (2009).
  * Trickey  _et al._ (2009) S. Trickey, V. Karasiev, and R. Jones, International Journal of Quantum Chemistry 109, 2943 – 2952 (2009), cited by: 25.
  * Huang and Carter (2010a) C. Huang and E. A. Carter, Phys. Rev. B 81, 045206 (2010a).
  * Karasiev and Trickey (2012) V. Karasiev and S. Trickey, Computer Physics Communications 183, 2519 (2012).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Parr and Yang (1994) R. G. Parr and W. Yang,  _Density-Functional Theory of Atoms and Molecules (International Series of Monographs on Chemistry)_ (Oxford University Press, USA, 1994).
  * Thomas (1927) L. H. Thomas, in  _Mathematical proceedings of the Cambridge philosophical society_ , Vol. 23 (Cambridge University Press, 1927) pp. 542–548.
  * Fermi (1927) E. Fermi, Rend. Accad. Naz. Lincei 6, 5 (1927).
  * Weizsäcker (1935) C. F. v. Weizsäcker, Journal of Physics 96, 431 (1935).
  * Teller (1962) E. Teller, Rev. Mod. Phys. 34, 627 (1962).
  * Wang and Teter (1992a) L.-W. Wang and M. P. Teter, Phys. Rev. B 45, 13196 (1992a).
  * Witt  _et al._ (2018a) W. C. Witt, B. G. del Rio, J. M. Dieterich, and E. A. Carter, Journal of Materials Research 33, 777–795 (2018a).
  * Wang  _et al._ (1999) Y. A. Wang, N. Govind, and E. A. Carter, Physical Review B 60, 16350 (1999).
  * Hung and Carter (2009) L. Hung and E. A. Carter, Chemical Physics Letters 475, 163 (2009).
  * Wang and Teter (1992b) L.-W. Wang and M. P. Teter, Phys. Rev. B 45, 13196 (1992b).
  * Huang and Carter (2010b) C. Huang and E. A. Carter, Phys. Rev. B 81, 045206 (2010b).
  * Mi  _et al._ (2018) W. Mi, A. Genova, and M. Pavanello, The Journal of Chemical Physics 148, 184107 (2018), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5023926/15540997/184107_1_online.pdf .
  * Shao  _et al._ (2021) X. Shao, W. Mi, and M. Pavanello, Phys. Rev. B 104, 045118 (2021).
  * Witt  _et al._ (2018b) W. C. Witt, B. G. del Rio, J. M. Dieterich, and E. A. Carter, Journal of Materials Research 33, 777–795 (2018b).
  * Snyder  _et al._ (2012) J. C. Snyder, M. Rupp, K. Hansen, K.-R. Müller, and K. Burke, Phys. Rev. Lett. 108, 253002 (2012).
  * Meyer  _et al._ (2020) R. Meyer, M. Weichselbaum, and A. W. Hauser, Journal of Chemical Theory and Computation 16, 5685 (2020).
  * Sun and Chen (2024) L. Sun and M. Chen, Phys. Rev. B 109, 115135 (2024).
  * Seino  _et al._ (2018) J. Seino, R. Kageyama, M. Fujinami, Y. Ikabata, and H. Nakai, The Journal of Chemical Physics 148, 241705 (2018), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5007230/16652841/241705_1_online.pdf .
  * Hollingsworth  _et al._ (2018) J. Hollingsworth, L. Li, T. E. Baker, and K. Burke, The Journal of Chemical Physics 148, 241743 (2018), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5025668/16656974/241743_1_online.pdf .
  * Tal and Bader (1978) Y. Tal and R. Bader, International Journal of Quantum Chemistry 14, 153 – 168 (1978), cited by: 153.
  * Bartolotti and Acharya (1982) L. Bartolotti and P. Acharya, The Journal of Chemical Physics 77, 4576 – 4585 (1982), cited by: 47.
  * Levy and Ou-Yang (1988) M. Levy and H. Ou-Yang, Physical Review A 38, 625 – 629 (1988), cited by: 180\.
  * Finzel (2021a) K. Finzel, Molecules 26, 1539 (2021a).
  * Finzel (2021b) K. Finzel, Acta Crystallographica Section B: Structural Science, Crystal Engineering and Materials 77 (2021b).
  * Alonso and Girifalco (1978) J. A. Alonso and L. A. Girifalco, Physical Review B 17, 3735 (1978).
  * Alonso and Balbás (1981) J. Alonso and L. Balbás, Physics Letters A 81, 467 (1981).
  * Deb and Ghosh (1983) B. M. Deb and S. K. Ghosh, International Journal of Quantum Chemistry 23, 1 (1983).
  * Deb and Ghosh (1982) B. M. Deb and S. K. Ghosh, The Journal of Chemical Physics 77, 342 (1982), https://pubs.aip.org/aip/jcp/article-pdf/77/1/342/18937204/342_1_online.pdf .
  * Deb and Chattaraj (1989) B. M. Deb and P. K. Chattaraj, Phys. Rev. A 39, 1696 (1989).
  * Roy  _et al._ (1999) A. K. Roy, B. K. Dey, and B. M. Deb, Chemical Physics Letters 308, 523 (1999).
  * Priya and Sadhukhan (2022) Priya and M. Sadhukhan, Molecular Physics 0, e2136114 (2022), https://doi.org/10.1080/00268976.2022.2136114 .
  * Economou (2006) E. Economou,  _Green’s Functions in Quantum Physics_, Springer Series in Solid-State Sciences (Springer, 2006).
  * Note (1) Consider Helium and Lithium atoms with exact non-relativistic energies as $\approx-2.904$ a.u. and $\approx-14.864$ a.u., respectively. To them the hydrogenic parts of the Hamiltonian attribute $\approx-4.0$ a.u. and $\approx-20.25$ a.u. The contributions of the inter-electronic repulsion are evidently smaller compared to the hydrogenic part.
  * Yang (1986) W. Yang, Physical Review A 34, 4575 (1986).
  * Brack  _et al._ (1976) M. Brack, B. Jennings, and Y. Chu, Physics Letters B 65, 1 (1976).
  * Dyson and Lenard (1967) F. J. Dyson and A. Lenard, Journal of Mathematical Physics 8, 423 (1967), https://pubs.aip.org/aip/jmp/article-pdf/8/3/423/19327819/423_1_online.pdf .
  * Lieb and Thirring (1975) E. H. Lieb and W. E. Thirring, Phys. Rev. Lett. 35, 687 (1975).
  * Casas  _et al._ (2012) F. Casas, A. Murua, and M. Nadinic, Computer Physics Communications 183, 2386 (2012).
  * Magnus (1954) W. Magnus, Commun. Pure Appl. Math. , 649 (1954).
  * Gell-Mann and Brueckner (1957) M. Gell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 (1957).
  * Feynman  _et al._ (2010) R. Feynman, A. Hibbs, and D. Styer,  _Quantum Mechanics and Path Integrals_, Dover Books on Physics (Dover Publications, 2010).
  * Chua  _et al._ (2018) C.-K. Chua, Y.-T. Liu, and G.-G. Wong, Journal of Physics Communications 2, 035007 (2018).
  * Acharya  _et al._ (1980) P. K. Acharya, L. J. Bartolotti, S. B. Sears, and R. G. Parr, Proceedings of the National Academy of Sciences of the United States of America 77, 6978 (1980).
  * Clementi and Roetti (1974) E. Clementi and C. Roetti, Atomic Data and Nuclear Data Tables 14, 177 (1974).
  * Kramida  _et al._ (2023) A. Kramida, Yu. Ralchenko, J. Reader, and and NIST ASD Team, NIST Atomic Spectra Database (ver. 5.11), [Online]. Available: https://physics.nist.gov/asd [2024, July 10]. National Institute of Standards and Technology, Gaithersburg, MD. (2023).
  * Wiles and Manby (2018) T. C. Wiles and F. R. Manby, Journal of Chemical Theory and Computation 14, 4590 (2018), pMID: 30080967, https://doi.org/10.1021/acs.jctc.8b00337 .
  * Lieb (1976) E. H. Lieb, Rev. Mod. Phys. 48, 553 (1976).