{'fraction_non_alphanumeric': 0.04735665644331237, 'fraction_numerical': 0.01174819358527837, 'mean_word_length': 4.466894003978403, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kinetic energy density for open-shell systems: Analysis and development of a novel technique

The quest for an approximate yet accurate kinetic energy density functional is central to the development of orbital-free density functional theory. While a recipe for closed-shell systems has been proposed earlier, we have shown that it cannot be na\"ively extended to open-shell atoms. In this present work, we investigated the efficacy of an ad-hoc recipe to compute the kinetic energy densities for open-shell atoms by extending the methodology used for closed-shell systems. We have also analyzed the spin-dependent features of Pauli potentials derived from two previously devised enhancement factors. Further, we have proposed an alternate but exact methodology to systematically compute the kinetic energy density for atoms of arbitrary spin multiplicity.


I. INTRODUCTION
Kinetic energy density (KED) functional is the toughest bottleneck for developing an effective orbital-free density functional theoretic framework for atoms and molecules. However, most of the developments in that direction have not considered the effects of spin-polarization explicitly.
In this work, we have extended some previously developed forms of kinetic energy density to open-shell systems and introduced a novel and exact route to derive kinetic energy density in a systematic manner.
The fundamental promise of density functional theory is that the total energy of an electronic ground state can be completely characterized by its one particle reduced density [1,2]. Despite this assurance, the analytical form of KED and inter-electronic repulsion potential as functionals of electron density is unknown. To bypass this problem, Kohn and Sham [3] proposed a scheme where the density is partitioned into orbitals and computed from the approximated inter-electronic repulsion and kinetic energy density. The problem of kinetic energy, though hidden in KS scheme, remains an open challenge to this day [4]. While in principle exact, Kohn-Sham equations cannot be applied to large systems due to prohibitive computational cost [5]. In OFDFT, the total electron density can be computed from an Euler equation [4,6,7] or imaginary time-evolution of a time-dependent hydrodynamical equation [8,9]. As a result, the computational cost does not depend upon the number of electrons in the system explicitly, making the search for an accurate OFDFT method important. Recent efforts to employ OFDFT methods with the machine learning algorithms [10,11] and pseudopotentials [12] method manifest the usefulness of an accurate and general form of KED for simulations of large scale materials.
Several attempts have been formulated over the years to obtain an accurate description of exact kinetic energy density [4]. For solid-state systems, the formulation of adequately accurate KED has been achieved [13,14]. However, an accurate KED for atomic and molecular systems which produces atomic shell structures via self-consistent field calculation is lacking [7,15]. Almost all approximations of KED are based on the Thomas-Fermi model [16,17]. For solid state systems generalized-gradient approximations (GGA) [18][19][20], meta-GGA [18,21], conjoint gradient correction approach [22], as well as non-local response-based approaches [23,24] enjoyed some success. Two points weighted density approximation of KED [25] and bi-functional-based methods are also being developed at present [26,27]. The other related quantity essential to self-consistently compute correct electron density is Pauli potential [28,29]. This quantity has long been associated with proper electronic shell structure [15,30,31]. Recently, shell-structure-based functionals are also showing promises [32,33].
We arranged the article as follows. First, we have introduced (1) an ad-hoc procedure, (2) a GGA-functional based procedure and (3) an exact methodology to compute the enhancement factors for spin polarized systems in section II. In section III we have rigorously re-parameterized and analyzed the Pauli potentials derived from GGA functionals for open and closed-shell systems.
The paper concludes with section IV.

II. THEORY
The kinetic energy [34] of an electronic system is given by where the KED is defined by (atomic units have been used throughout the article) For a single determinantal wave function the one particle density matrix can be written as Here ε α and φ α i are highest occupied level of α spin orbitals and i th orbitals of α spin manifold, respectively. If we use Eq.(4) for density matrix, where the orbitals are Kohn-Sham orbitals, then we get the corresponding kinetic energy T s [ρ].
In orbital-free density functional theory (OFDFT), the total kinetic energy functional is written as, Here, the first term is the von Weizsäcker kinetic energy. The second term is known as modified Thomas The Hartree-Fock KED t HF is computed using the Hartree-Fock orbitals in Eq.(4).
Note that, the Kohn-Sham limit of kinetic energy is not the true kinetic energy and a part of it is dumped in exchange-correlation functional. In terms of Levy-Lieb constrained-search formulation of DFT, we can construct the Kohn-Sham equation from defining the extremum condition for T s .
This approach, however, gives rise to Wang paradox [4]. Nevertheless, in this article, we will use Pauli potential as a shorthand for T corr .
For self-consistent calculation of electron density, one can use the imaginary-time evolution of Here the hydrodynamical function ψ(r,t)yields the electron density as ρ(r) = |ψ(r)| 2 . V Coul (r) and V XC (r) are electrostatic (of classical mechanical origin) and exchange-correlation (of quantum mechanical origin) potentials, respectively. The last term in Eq.(9) is known as the Pauli potential where C T F = 3 10 (3π 2 ) 2/3 . In their works, Deb et. al.. modeled g(r) as a sum of Gaussian functions.
They have used one Gaussian function per shell. The same technique has recently been used by Finzel [7] to model the Pauli potential where the Gaussian functions were replaced by Slater-type functions. Recently efforts have been made to approximate the Pauli potential via a Padé-type functional of r as well as ρ(r) [32].
Eq.(5) and subsequent derivation of Eq.(8), however, assumes that the two-particle density matrix can be written as The formulation of F(r) is then followed by defining the correlation factor such that one can compute one-particle density matrix as For homogeneous electron gas, the correlation factor is exactly known. A modification due to inhomogeneous electron density results a modified corrlation function, which in turn yields F(r).
However, Eq.(12) and Eq.(11) are strictly correct for closed-shell systems. For an open shell system the two-particle density matrix for a single determinantal wave function can be written as [6] ρ 2 (r, Eq.(15), unlike Eq.(11), cannot be inverted to obtain the expressions for ρ α 1 (r, r ′ ) or ρ β 1 (r, r ′ ) separately. Moreover, since due to Eq.(3), we cannot use Eq.(13) to define the correlation function for open-shell systems.
Due to the same reasons, a straightforward extension of spin-polarized Pauli potential is also unavailable.
A. An ad-hoc definition of spin-polarized enhancement factors Extensions of Eq.(7) and Eq. (10) for open-shell systems, are therefore, necessary for a generalized description of atoms and molecules. The first attempt, albeit an ad-hoc one, is inspired by the spin-polarized extension of von Weizsäcker and Thomas-Fermi kinetic energy [6] given by and respectively. We therefore define spin-polarized Pauli kinetic energy where spin variable σ ∈ {α, β }. Note that the factor 2 2/3 is absorbed inside F σ (r). The spinpolarized Pauli potential is therefore can be computed by B. Generalized gradient approximation-based Pauli potentials Different classes of enhancement factors [22] have been introduced over the years. In a recent work [14] some of the enhancement factors based on generalized gradient approximation (GGA) have been analyzed for semiconductors and metals. Here we will consider two of their simplest enhancement factors and analyze the corresponding Pauli potentials for open and closed-shell systems. These two enhancement factors are given by

the linear function
and 2. the Pauli-Gaussian(PG) enhancement function where λ and µ are free parameters. Both of them have been expressed in terms of scaled gradient of density where C 0 = 2(3π 2 ) 1/3 . Other forms of F(r) and further modifications by introducing scaled Lapla- have also been considered recently [21,36,37]. However, we will not consider them in this work.
The Pauli potential and therefore the corresponding g(r) can be computed by The Pauli potentials corresponding to Eq. (22) and Eq. (23) are where p = s 2 . The µ in Eq. (27) is defined as µ = 5 3 (λ − 1). Note that s is a dimensionless quantity. These expressions can therefore be extended to spin-polarized cases phenomenologically by replacing ρ(r) by ρ σ (r) and s by s σ 2 −1/3 . Here C. An exact and systematic method using Green's function Now we will lay an exact road map for deriving the enhancement factor and therefore Pauli potential using Green's function (GF) method. We start from defining the Green's function G(r, r ′ ; β ) as the Laplace transform of the density matrix [6] G(r, Since the Laplace transform is a linear transform, we can extend the definition of Green's function to spin-polarized cases in a straightforward manner (cf.Eq. (3)). We can obtain density matrix for a given spin multiplicity via Bromwich integral where γ ∈ R > . As a result, G σ (r, r ′ ; β ) can be written as whereĤ is the corresponding Hamiltonian with eigenvalues {ε σ i } and eigenfunctions {φ σ i }. We will drop the superscript σ for Green's functions in this section for brevity. However, every Green's function discussed onward will be assumed as spin-polarized.
We can write the Hamiltonian of an atom aŝ whereĤ 0 ,Ĥ Z andŴ are free-particle Hamiltonian, electron nucleus attraction potential and inter-electronic repulsion potentials respectively.Ĥ H =Ĥ 0 +Ĥ Z is the sum of Hamiltonians of a Hydrogenic atom with Z. Consequently, the density matrix corresponding toĤ H produces von Weizsäcker kinetic energy with correction which can be computed in a systematic manner [38].
We define two more Green's functions and corresponding to the Hydrogenic system and free-particle systems respectively. Note thatĤ HEG = H 0 +Ŵ describes the homogeneous electron gas limit and we will call the corresponding Green's function G HEG (r, r ′ ; β ).
A connection between these Green's functions can be obtained Using Zassenhaus formula [39,40] we can writê We find since Ĥ Z ,Ŵ commutes. As a result, all commutator in Eq. (38) can be replaced by Ĥ 0 ,Ŵ .
As a result, Using Eq. (40) in Eq. (36) and employing resolution of identity, we obtain Eq.(41), similar to a Dyson equation, presents the total Green's function as an infinite series and is an exact result, leading to a systematic way to compute the total Green's function. Using the definition of Green's functions and considering only leading order term (n = 0) of Eq. (41), we obtain Since bothĤ H andĤ HEG can be written asĤ 0 + other terms , we can again apply Zassenhauss formula for G H (r, r ′ ; β ) and G HEG (r, r ′ ; β ) to obtain Green's functions as G HEG (r, r ′ ; β ) = G 0 (r, r ′ ; β ) + . . .
Using Eqs. (43) and (44) in Eq. (42), we obtain Using the definition of Green's function used here we find therefore the leading order term of the Green's function is It can be seen from Eq.(31) that Again using the leading-order term of Eq.(47) in Eq.(46), we find that where ρ H 1 (r, r ′ ) and ρ 0 1 (r, r ′ ) are density matrices for a Hydrogenic atom and a free-particle system, respectively. Clearly, the form of kinetic energy from this approach turns out to be The higher order terms now can be considered systematically. However, such a task will require significant mathematical exercise and will be considered later. Eq.(49) is the form described by Carter et. al. [13,31].

III. RESULTS AND DISCUSSION
We have considered two noble gas atoms (Neon, Argon) and four first-row elements (Lithium, Beryllium, Boron and Carbon) as our representative systems. Three of them (Ar, Ne, Be) are spin-paired while the other three (Li, B, C) are spin-polarized. Among them, Li also has zero beta spin in second shell. This variety of electron structure is necessary to understand the effects of spin-multiplicity on the Pauli potential. For larger atoms, more intricate effects such as spin-orbit coupling and relativistic corrections become important as well. These effects are outside of the scope of this work. We will consider them in future works.
All densities are computed from parametrized Slater functions for the Hartree-Fock wave functions [41]. All benchmark results for F(r) and g(r) (dubbed F HF (r) and g HF (r), respectively onward) are computed using the parameters provided in Ref. [9].
A. Optimization of free parameters for g E (r) and g L (r) We begin by optimizing the value of parameter µ used in Eqs. (27) and (28). To accomplish that, we have computed a metric of deviation from the benchmark g HF (r) as A value of µ which minimizes E(µ) can be chosen as the optimum choice. We have chosen Ne and Ar atoms as our benchmark systems. Accurate g HF (r) for them are available as a linear combination of parametrized Gaussian functions [9,30]. These parameters were set to reproduce the Hartree-Fock level electron density upon self-consistent field calculations.
Variations of E(µ) as a function of µ are presented in Fig.1.
For Ar atom, the minima are clear for both g L (r) and g E (r). For Ne, such clear minimum is not pronounced. Therefore, we have chosen the value for µ which minimizes E(µ) for Ar atom. It is clear that the parameter µ is actually system-dependent and therefore is not transferable. A better understanding of system dependence of µ is required if these enhancement factors are to be used for atomic systems. For all our subsequent results, µ = 2.5 for g L (r) and µ = 0.75 for g E (r) have been used irrespective of the systems.

B. Comparison between different g(r) for closed-shell systems
Having optimized the parameters, we now examine the features of each g(r) for closed shell atoms. We have compared them with the benchmark g HF (r) functions mentioned above (Fig.2).
In this plot both g L (r) and g E (r) do not behave well for r → 0. A possible reason for this behavior is the Laplacian terms in Eq. (27) and Eq. (28). g L (r) (black solid line) mimics the shell structure for both Ne and Ar atoms more accurately. However, for large r, the function diverges quickly.
Here large values of r signify the distances where both the electron density and its gradients are very small. As a result, our computed g(r) is also susceptible to numerical inaccuracies. On the other hand, g E (r) (blue solid line) shows very different profile while it is more well-behaved at large values of r. It appears that, although simpler in form, g L (r) produces better shell-structures than the more complicated g E (r) contrary to bulk systems [14,42]. for open-shell atoms. For such systems, we examined g L (r), g α L (r) and g β L (r) for the total density, α spin density and β spin density, respectively (Fig.3). For all four cases, g β L (r) showed the highest peak followed by g α L (r) and g L (r). This trend indicates that the peak height is inversely proportional to the number of electrons considered. This is also corroborated by Fig.2 which clearly shows that Neon has the lowest peak height among all first-row elements considered. The peak height corresponding to the inner-shell of Argon atom is even lower supporting our observation.
Also, the peak positions shift to smaller distances as the number of electrons increases, showing the contraction of electron density with increasing nuclear charge [? ]. As mentioned before g L (r) diverges for large value of r, explaining the diverging g β L (r) (no β electron in 2s) for Li. For Be the g α L = g β L , as expected. B and C have same number of β electrons and hence exhibit very similar g β L . For C atom, g α L (r) is very similar to that of g L (r). While the shell structure is clearly visible for all four atoms, note that the g L values sometimes fall below zero, rendering the corresponding Pauli potential negative. This is a serious drawback of this g L (r). Next we examine the variation of peak heights of g L (r) with atomic number Z (same as electron number N here) for first-row elements (Fig.4). As expected, the peak height decreases with Z. To Next, we analyze the effects of g E (r) for open-shell atoms (Fig.5). Similar to the cases discussed in the previous section, we considered g E (r), g α E (r) and g density and β spin density, respectively. Following the same trend already observed in sectionIII B, we find that g E (r) is well-behaved for large r. However, the absence of any clear shell structures, akin to g L (r), obscures the applicability of Eq. (28) in Pauli potential. Unlike the divergence observed for g β L (r) for Li, here we see that the g β E (r) settles to zero for r ≥ 1.2. The problem of negative values for g E (r) is more prominent here compared to g L (r). The other features are shown by g α E and g β E are not very clear at present. A more careful analysis of Eqs. (27) and (28) is therefore required to develop better Pauli potentials in the future.

IV. CONCLUSION
In summary, we have (1) pointed out the mathematical difficulty to extend a previously devised derivation for enhancement factors for open-shell atoms, (2) analyzed Pauli potentials derived from two previously devised enhancement factors, and (3) developed an exact and novel method to compute the enhancement factors which do not depend on the spin-multiplicity of the system.
To achieve this we employed Green's function technique to obtain a well-defined but infinite series expression for density matrix. We found that the Pauli potentials derived from previously used enhancement factors does not always meet the expected asymptotic behaviors. Also, the linear enhancement factor produces the atomic shell structure better compared to the more complicated Gaussian enhancement factors. However, both of these functions fail to provide a correct descrip- tion of g(r) near the nucleus. Further studies are required to understand the reasons and remedies for this problem. Moreover, we have found that the shell structure in g(r) corresponding to linear enhancement factor is very sensitive to spin multiplicity. This work raises some new questions as well. The possibility to derive adequately accurate yet analytically closed forms of enhancement factor and Pauli potential is the most significant among them. While the Green's function method is promising, the convergence of the infinite series remains an open question. Description of Pauli potential for molecular systems is another problem that may be taken up in the future. Furthermore, the effect of the bond-breaking and bond-making process on these quantities are yet to be explored. We will be addressing these questions in our future works.