Path Integrals for Electronic Densities, Reactivity Indices, and Localization Functions in Quantum Systems

The density matrix theory, the ancestor of density functional theory, provides the immediate framework for Path Integral (PI) development, allowing the canonical density be extended for the many-electronic systems through the density functional closure relationship. Yet, the use of path integral formalism for electronic density prescription presents several advantages: assures the inner quantum mechanical description of the system by parameterized paths; averages the quantum fluctuations; behaves as the propagator for time-space evolution of quantum information; resembles Schrödinger equation; allows quantum statistical description of the system through partition function computing. In this framework, four levels of path integral formalism were presented: the Feynman quantum mechanical, the semiclassical, the Feynman-Kleinert effective classical, and the Fokker-Planck non-equilibrium ones. In each case the density matrix or/and the canonical density were rigorously defined and presented. The practical specializations for quantum free and harmonic motions, for statistical high and low temperature limits, the smearing justification for the Bohr’s quantum stability postulate with the paradigmatic Hydrogen atomic excursion, along the quantum chemical calculation of semiclassical electronegativity and hardness, of chemical action and Mulliken electronegativity, as well as by the Markovian generalizations of Becke-Edgecombe electronic focalization functions – all advocate for the reliability of assuming PI formalism of quantum mechanics as a versatile one, suited for analytically and/or computationally modeling of a variety of fundamental physical and chemical reactivity concepts characterizing the (density driving) many-electronic systems.


Introduction
In modern conceptual and computational chemistry the Density Functional Theory (DFT) [1][2][3][4][5][6][7][8][9][10][11][12] plays the central role since its capabilities in providing both structural and reactivity information about the atoms and of their interaction in molecules and nanostructures [13][14][15][16][17][18][19][20]. Yet, the main vehicle stays the electronic density, an observable quantity, which intimately relates with the more abstract quantum mechanical concept of wave function through the basic relationship [3] written for a collection of N-many-electronic system with space-spin coordinates Besides the sub-script index, it is worth noting that with Equations (1) and (2) either the ground state or the valence state(s) of a many body system may be computed by applying the variational principle upon the total or valence density functional energy [7] [ , and of the so called chemical action term [7][8][9][10] [ ] ( ) ( )dx accounts in principle for all observable effects an electronic structure may manifest as the chemical reactivity. Yet, although accessible experimentally [21], for a deeper comprehension of the physical-chemical phenomenology of bonding the electronic density should be "visualized" also within an analytical framework. In this regard the already consecrated attempt of density matrix made history over a half century by employing the so called ensemble (statistical) density operator [22][23][24][25][26]   ; x x ρ provides by integration the total number of electrons in the same way as the observable density does [27][28][29][30][31] ( ) that is formally equivalent with the "trace" operation on the bilocal density matrix ( )  (8) From this perspective it is clear that having in hand a reliable method to express, for various applied potentials ) (x V the bilocal density matrix, from here on called simple as density matrix, yields in fact the electronic density itself. Fortunately the quantum mechanically formalism had advanced a rigorous way of expressing the density matrix since recognizing it as a specialization of the more general time evolution quantum amplitude, a.k.a. the quantum propagator or the Green function [32][33][34][35][36] ( ) Moreover, within the thermal-temporal (quantum mechanical-to quantum statistical) Wick transformation ; exp exp (11) to release the uni-particle density while fulfilling the normalization constraint given by Equation (2). The competition between the variation of density functional of energy (3) and the density itself (12), i.e., the ratio between the global and local influences on a many-electronic system (equilibrium, aromaticity, etc.) is quantified by the modern electronegativity index and by its chemical hardness companion [3,[5][6][7][8][9][10][11][12], given respectively as which may assume various analytical realizations within density functional theory [7,9], while playing the crucial roles in driving chemical bonding and reactivity [12,20]; this may be also immediately seen from Equations (13) when electronegativity is further identified with the negative of the chemical potential of a given system ( ) [6]: if electronegativity is assimilated with the chemical potential, the chemical hardness -as its derivative -corresponds with the chemical force; thus, they together constitute the minimal necessary set of indices to consistently model the complete scenario of chemical reactivity, from encountering adducts to the stabilized products [9]. For these reasons they will be systematically presented in this review as application for different levels of path integral approximations in either matrix density or density computations. However, it is worth comment that the above factorization only apparently assumes the manyparticle system without internal interaction (exchange and correlation), while all these effects are to be incorporated in the way the bare applied potential is replaced with an effective one or by performing a variational (upon the) perturbation procedure for optimizing the (bilocal) density matrix in accordance with above recipe. Equation (14) is nevertheless nothing else than the so called path integral of the density matrix or of the time evolution quantum amplitude, while ∫ ) (τ Dx is a complex symbol of integration over (parameterized, quantum) paths that when unfolded recovers the normal shape of integration, i.e., with separated variables of integration, e.g., Equations (117)

or (316) and
Refs. [37][38][39]. The way in which this path integral is computed for various working potentials will generate the analytical solution of the quantum amplitude, and implicitly to the density matrix, from where the electronic density is immediately found by the identity ( ) ( ) . Then, having the electronic density any known or approximated density functional may be evaluated and employed in describing the chemical structure and reactivity in an analytical manner, while allowing better conceptual understanding of the obtained models and predictions. With these, the motivation behind the present project becomes clear: computing the electronic density relays in fact on evaluating the associate path integral for a given potential; moreover, the many-body effects are to be resumed in rewriting the applied potential into an effective one (one way) or to advance the so called variational perturbation methodology in order the inter-particle effects be accommodated [39][40][41][42][43][44][45][46][47][48], being this latter approach left for further communications.
Consequently, the review unfolds on bigger scale the ideas here presented: it starts with the basic properties of the density matrix and showing how the path integral concept arises naturally in this framework. Then, a more formal introduction of path integral methodology is presented in the spirit of Richard Feynman, its main promoter; and the use of path integrals is exemplified in computing semiclassical time evolution amplitudes with application on atomic electronegativity and chemical hardness reactivity indices. The simplified many-body approach is then given through exposing the Feynman-Kleinert algorithm for effective potentials, with application on computing atomic Mulliken electronegativities, while the non-equilibrium Fokker-Planck approach is exposed and applied in the context of Markovian stochastic motion within the anharmonic potential and then extended to modeling the electronic localization through computing several Markovian electronic functions while comparing them with the circulating Becke-Edgecombe one [49,50]. This way, a fruitful step is hopefully made towards unifying the physical-chemical principles of electronic structure and reactivity on a meaningful quantum basis.  (21) with the sum of diagonal matrix elements (yielding the "trace" function)  (23) Note that through the above deductions the double (independent) averages technique was adopted exploiting therefore the associate sums inter-conversions to produce the simplified results [51][52][53]. Yet, this technique is equivalent with quantum mechanically factorization of the entire Hilbert space into sub-spaces or, at the limit, into the subspace of interest (that selected to be measured, for instance) and the rest of the space, being thus this approach equivalent with a system-bath sample; this is useful for better understanding the stochastic phenomena -to be latter exposed -that underlay to open quantum systems. Therefore, such mechanism may be considered as belonging to the physical foundations for the chemical reactivity.

On Mono-, Many-, and Reduced-Electronic Density Matrices
Next, in the case the concerned quantum states are eigen-states, they fulfill the normalization constraint  (24) on which base the above density operator now reads with the same form as presented din Introduction, Equation (6), from where there appears that the eigen-equation for it looks like giving with the eigen-values (as the diagonal elements) as k k k w = ϕ ρ ϕˆ (26) as the observed values of the averaged density operator. Since they are weights of probability they have to naturally fulfill the closure probability relationship over the entire sample 1 = ∑ k k w (27) from where follows the "normalization of density operator" through its above Trace property of Equation (22) 1 = ρ Tr (28) Moreover, in these eigen-conditions, the operatorial average further reads from Equation (23) ( ) Now, there appears with better clarity the major role the density operator plays in quantum measurements, since it convolutes with a given operator to produce its (averaged) measured value on the prepared eigen-states. Nevertheless, when the so called pure states are employed or prepared, the preceding distinction between the subsystem and system vanishes, and the density operator takes the pure quantum mechanical form of an elementary projector Λ ≡ =φ ϕ ρ (30) This is a very useful expression for considering it associated with the mono-density operators when many-fermionic systems are treated, although a similar procedure applies for mixed (sample) states as well. It is immediate to see that for N formally independent partitions the Hilbert space corresponding to the N-mono-particle densities on pure states, we have through Equations (22), (27), (28) and (30) producing the total operator -projector constructed by the sum while correctly normalized to the total number of particles r Tr Tr (33) Yet, the anti-symmetric restriction the N-fermionic state may be accounted from the monoelectronic states through considering Slater permutated ( α P ) products [8,35] (34) for constructing the N-electronic density operator  (36) However, in practice, due to the fact the multi-particle operators associate with number of systemic properties less than the total number of particle, say of order N p < , worth working with the p-order reduced density matrix introduced as ρ (37) with the following features [27].
where the first order density matrix casts as abstracted from general definition With these concepts it is worth noting the major importance that the first order density plays in computing the higher order reduced density matrices that in turn enter the operatorial averages, for instance A special reference may be made in regard of the free-relativist treatment of many-electronic atoms, ions, bi-or poly-atomic molecules governed by the working Hamiltonian  (43) whose terms are represented the inter-nuclear repulsion (only for molecules), free electronic motion, electron-nuclei Coulombic attraction, and inter-electronic Coulombian repulsion, respectively. For such Hamiltonian the average value is computed through considering electronic density of the first or second order only there where the electronic influence is present, while the degree of matrix density is fixed by the type of electronic interaction It is obvious that although the second order reduced matrix has appeared, its general form ρ (45) may be further reduced to the first one through the above determinant rule, see Equation (40) this way emphasizing on the importance of the first order reduced matrix knowledge. The astonishing physical meaning behind this formalism relays in the fact that any multi-particle interaction (two-particle interaction included) may be reduced to the single particle behavior; in other terms, vice-versa, the appropriate perturbation (including strong-coupling) of the single particle evolution carries the equivalent information as that characterizing the whole many-body system.
In fact, the power of the density matrix formalism resides in reducing a many-body problem to the single particle density matrix, abstracted from the single Slater determinant of Equation (36) (48) considerably simplifies the quantum problem to be solved. Let's illustrate this by firstly quoting that Fock-Dirac density operator of Equation (48) has two fundamental properties, namely: o The idempotency Remarkably, the last two identities may serve as the constraints when minimizing the above Hamiltonian average, here appropriately rewritten employing Equations (44) and (46), and where all external applied potential, were resumed under generic ) ( 1 x V quantity producing the actual so called Performing the functional derivative respecting the Fock-Dirac electron density in (54) (55) which eventually transcribes at the operatorial level as meaning that they both admit the same set of eigen-functions. This is nevertheless the gate for obtaining the density (matrix) functional energy expressions by means of finding the density (matrix) eigen-solutions only.
Yet, condition (59b) is indeed a workable (reduced) condition raised from optimization of the averaged Hamiltonian of a many-electronic system, since the more general one referring to the whole Hamiltonian, known as the Liouville or Neumann equation, is obtained employing the temporal Schrödinger equation: to the evolution equation of Fock-Dirac density operator evolution Lastly, it should be noted that all above properties may be rewritten since considering the mixed porder reduced matrix with the form as a natural extension of the pure states. However, the sample statistical effects may be better considered by further expressing the electronic density operator and its matrix, the associate equation and the properties for systems in thermodynamic equilibrium (with environment) -a matter addressed in next section.

Canonical Density, Bloch Equation, and the Need of Path Integral
For a quantum system obeying the N-mono-electronic eigen-equations the probability of finding one particle in the state k ϕ at thermodynamical equilibrium with others, while all states are considered as a closed supra-system with no mass or charge transfer allowed, is given by the canonical distribution [36] ( ) ( ) (63) providing the mixed Fock-Dirac density with the form This is a very interesting and important result motivating the quantum statistical approach in determining the density of states since it corresponds to the N-sample particles throughout a simple Nmultiplication. Note that Equation (64) is in full agreement with that introduced in Equation (12), and very well suited for handling since respecting the DFT custom normalization of Equation (2), while its normalization factor, the partition function ( ) , follows from such constraint with the consecrated The recognized importance of partition functions in computing the internal energy as the average of the Hamiltonian of the system: Tr Tr (66) or to evaluate the free energy of the system is thus transferred to the knowledge of the closed evolution amplitude x e x Ĥ β − , that at its turn is based on the genuine (not-normalized) density operator sometimes called also like canonic density operator.
The great importance of density operator of Equation (68) is immediately visualized in three ways: o It identifies the evolution operator on the ground of Wick equivalence relationship of Equation (10), which allows the transformation of the Schrödinger into Heisenberg or Interaction pictures for appropriately describing the quantum interactions [53]; o It produces the so called Bloch equation [21] by taking its β derivative that identifies with the Schrödinger equation for genuine density operator through the same Wick transformation given by Equation (10), thus providing the quantummechanically to quantum-statistical equivalence; o Fulfills the (short times, higher temperature) so called Markovian limiting condition a very useful constraint for developing either the perturbation or the variational formalism respecting electronic density and/or partition function, see below.
In the frame of coordinate representation the Bloch problem, i.e., the differential equation together with the initial (Cauchy) condition, looks like Solution of this system is a great task in general, unless the perturbation method is undertaken for writing the Hamiltonian as the sum of free and small interaction components for which the free Hamiltonian solution is completely known, say In these conditions, one may firstly write  (75) where the inter-Hamiltonian components were considered to freely commute as per wish; then, the that may be rearranged under the perturbative fashion correspondingly written in coordinate representation Such slicing procedure in solving the Bloch equation (72) for canonic density solution (80) seems an elegant way of avoiding the self-consistent equation (77). Therefore, it may further employed through reconsidering the problem (72) in a slightly modified variant, namely within the temporal approach that being of exponential type allows for direct slicing through factorization. That is, when considering the space partition given by coordinate cuts of (81), and assuming that the times flows equally on each sub-interval in quota of ε , , the density solution (83) may be written as a product of intermediary solutions towards the path integral representation (84) where the chained covariant density product was introduced along the extended integration metric The general canonic solution (84) is viewed as the path integral solution for the Bloch equation (82), being therefore as a necessity when looking to general solutions for a given Hamiltonian; it gives the general solution for electronic density (68) since accounting for all path connecting two end-points either in space and time (or temperatures) through in principle an infinite intermediary points; this way the resulted path integral comprises all quantum information contained by the particle' evolution between two states in thermodynamical equilibrium with environment (or the other mono-particle states). However, once having the canonical density evaluated by its path integral, the associate mixed density matrix may be immediately written employing the operatorial form (64) to the actual spatial representation with the path integral based partition function written in accordance with Equation (65) assuring the preservation of the general DFT normalization condition This way, the general algorithm linking the path integral to the density matrix and to the electronic density, most celebrated DFT quantity in computing various density functionals (energies, reactivity indices) for characterizing chemical structure and reactivity -was established, while emphasizing the basic role the path integral evaluation has towards a conceptual understanding of many-electronic quantum systems in their dynamics and interaction.
Being thus established the role and usefulness of path integral in density functional theory the next section will give more insight in appropriately defining (constructing) path integral such that to further facilitate its practical evaluation for electronic systems of physical-chemical interest.

Construction of the General Path Integral
Through reconsidering the slicing of (81) also for the time with the spatial ending points recalled as the quantum propagator of Equation (9), within the Wick equivalence (10) may be firstly rewritten in terms of associate evolution operator to successively become (93) where the last form was obtained when n-times the complete eigen-coordinate set was introduced for each pair of events with the elementary propagator between them ( ) ( ) on the elementary time interval Now, the elementary quantum evolution amplitude (95) is to be evaluated, firstly by reconsidering the eigen-coordinate unitary operator, in the working form to separate the operatorial Hamiltonian contributions to the kinetic and potential ones yielding: where we have used the first order limitation of the Baker-Hausdorff formula (100) by assuming the second order of elementary time intervals as vanishing Next, each obtained working energetic contribution is separately evaluated: for kinetic contribution the insertion of the momentum complete eigen-set while for potential elementary amplitude one gets With relations (103) and (104) Replacing the elementary quantum amplitude (105) back into the global one given by Equation (93), it takes the form which for an infinitesimal temporal partition, i.e. 0 ; → ∞ → ε n (107) behaves like the Feynman path integral of the quantum propagator [37,40,41,[54][55][56][57][58][59][60] ( ) through considering the limiting notations for the path integral measure and for the involved action Note that the results (108)- (110) give similar quantum information for the quantum evolution of a system as previously found with the mean of density matrix (84)- (86), yet in a more formal and general way throughout accounting all histories (possibilities for linking two events in time-space) for a quantum evolution [61][62][63] ( ) thus being suitable to be implemented in the N-particle density functional scheme (87)-(89) once it is analytically computed. For achieving such goal, a more practical form of the Feynman integral may be obtained once the Hamiltonian is implemented as leaving the action (110) unfolded as from where the momentum integrals in (109) is immediately solved to be (114) by formally applying the Poisson formula The remaining quantum evolution amplitude reads as the spatial path integral only assuming the actual modified measure of integration (117) and the working action Note that when the partition function (88) is under consideration, other path integral out of (116) has to be introduced by means of closed space coordinates, namely noting the new integration measure Therefore, at the first instance, some of the main advantages dealing with path integrals relay on following features: o Attractive conceptual representation of dynamical quantum processes without operatorial excursion; o Allows for quantum fluctuation description in analogy with thermic description, through changing the temporal intervals with the thermodynamical temperature by means of Wick transformation (10), i.e., transforming quantum mechanical (QM) into quantum statistical (QS) propagators from where the immediate writing of the associate QS-partition function having both QS objects written as the effect of transforming the canonical Lagrangean of action into the so called Euclidian one analogously with the fact the Euclidian metric has all its diagonal terms positively defined. Yet, the connection of the path integrals of propagators with the Schrödinger quantum formalism is to be revealed next.

Schrödinger Equation from Path Integral
There are two ways for showing the propagator path integral links with Schrödinger equation.

Propagator's Equation
Firstly, by employing one of the above path integral, say that of Equation (116) with (118) to perform the derivative Similarly for the second derivative we have while for time derivative we obtain by recalling the Hamilton-Jacobi equation of motion in the form Now, there is immediate that for a Hamiltonian of the form (112) one gets through multiplying both its side with the propagator (122) and then considering the relations (124) and (126), respectively, one leaves with the Schrödinger type equation for the path integral Remarkably, besides establishing the link with the Schrödinger picture, equation (127) tells something more important, namely that the wave function itself, i.e., ( ) , may be replaced (and generalized as well) by the quantum propagator ( ) ; this has a crucial consequence since the propagator is providing the N-electronic density in the directly and elegantly manner prescribed by the algorithm (87)- (89), here actualized as with partition function given as in (119), assuring for the correct N-representability (DFT) constraint: thus nicely replacing the complicated many-body wave function calculations.
Nevertheless, the path integral formalism is able to provide also the exact Schrödinger equation for the wave function, as will be shown in the sequel.

Wave Function's Equation
The starting point is the manifested equivalence between the path integral propagator and the Green function, with the role in transforming the wave-function registered on a space-time event into other one, either in the future of past quantum evolution. Here we consider only retarded phenomena modeled by the propagator ( ) ( ) which in accordance with the very beginning path integral construction, the slicing (90) and the relation (91), implies the existence of the so called quantum Huygens principle of wave-packet propagation [64] ( ) ( ) ( ) Yet, we will employ Equation (131) for an elementary propagator, modeling the quantum evolution presented in Figure 1, thus behaving like where A plays the role of the normalization constant in (132a) to assure the convergence of the wave function wave-packet. Equation (132a) may be still transformed employing the geometrical relation to compute the space and velocity averages 2 2 respectively, while changing the variable to furnish the actual form where Lagrangean was considered with its canonical form, as in (118), and the new constant factor was considered assimilating the minus sign of (133d). Next, since noting the square dependence of ξ in (132b) there will be assumed the series expansion in coordinate (ξ ) and time ( ε ) elementary steps restrained to the second and first order, respectively, being the time interval cut-off in accordance with the general (101) prescription. Thus we firstly have and the form (132b) successively rearranges: where the mixed orders producing a total order beyond the maximum equal two have been neglected, e.g., 0 2 ≅ εξ , and were we arranged the exponentials under integrals of Gaussian type (i.e., employing the identity ). Now, the integrals appearing on (135) are of Poisson type of various orders and, assuming the notation: are solved as: With these the expression (135) simplifies to which in the limit 0 → ε , commonly for path integrals, leaves with identity from where the convergence constant of path integral (132b) is found recovering the previous form, see Equation (114), thus confirming the consistency of the present approach. Nevertheless, with the constant (139b) back in (138) we get the equivalent forms being the last one identical with the Schrödinger wave function equation.
Thus it was therefore thoroughly proven that the Feynman path integral may be reduced to the quantum wave-packet motion while carrying also the information that connects coupled events across the paths' evolution, being by all of these a general approach of quantum mechanics and statistics.
The next section will deal with presenting practical application/calculation of the path integrals for fundamental quantum systems, e.g., the free and harmonic oscillator motions.

Path Integrals' Properties
There are three fundamental properties most useful for path integral calculations [65]. I. Firstly, one may combine the two above Schrödinger type bits of information about path integrals: the fact that propagator itself ( ) obeys the Schrödinger equation, see Equation (127), thus behaving like a sort of wave-function, and the fact that Schrödinger equation of the wavefunction is recovered by the quantum Huygens principle of wave-packet propagation, see Equation (131). Thus it makes sense to rewrite Equation (131) with the propagator instead of wave-function obtaining the so called group property for propagators which, nevertheless, may be recursively applied until covering the entire time slicing of the interval as given in (90) while remarking the absence of time intermediate integration.
II. Secondly, from the Huygens principle (131) there is abstracted also the limiting delta Diracfunction for a propagator connecting two space events simultaneously ( ) ( ) that is immediately proofed out This property is often used as the analytical check once a path integral propagator is calculated for a given system.

III.
Thirdly, and perhaps most practically, one would like to be able to solve the path integrals, say with canonical Lagrangean form (121a), in more direct way than to consider all multiple integrals involved by the measure (117).
Hopefully, this is possible working out the quantum fluctuations along the classical path connecting two space-time events. In other words, worth to disturb the classical path ) (t x cl by the quantum fluctuations ) (t x δ to obtain the quantum evolution path Very important, note that the quantum fluctuation vanishes at the end-points of the evolution path since "meeting" with the classical (observed) path, see Figure 2 ) being these constraints known as the Dirichlet boundary conditions. Now, one aims to separate the classical by the quantum fluctuation contributions also in the path integral propagator. Fortunately, this is possible for enough large class of potentials, more precisely for quadratic Lagrangeans of general type ( ) Actually, expanding the path integral action (118) around the classical path requires the expansion of its associate Lagrangean (146); so we get accordingly With the action (147b) one observes it practically separates into the classical and quantum fluctuation contributions; this has two major consequences: o Looking at the terms appearing in the whole Lagrangean (146) and to those present on the factor (150) it seems that once the last is known for a given Lagrangean, say L , then the same is characterizing also the modified one with the terms that are not present in the forms (150), namely o The resulting working path integral of the propagator now simply reads and gives intuitive inside of what path integral formalism of quantum mechanics really does: corrects the classical paths by the quantum fluctuations resumed as the amplitude of the (semi) classical wave. Next, the big challenge is to compute the above fluctuation factor (150); here there are two possible approaches. One is considering the fluctuations as a Fourier series expansion so that directly (although through enough involving procedure) solving the multiple integrals appearing in (150). This route was originally proposed by Feynman in his quantum mechanically devoted monograph [41], and recently refined by Kleinert in an extended textbook [39].
The second way is trickier, although with limitations, but it avoids performing the direct integration prescribed by (150), while being instructive since computing the quantum fluctuation again in terms of classical path action [65], however through employing the present first two propagator properties, the group property (141) and the delta-Dirac limit (143), upon the quantum wave (152).
As such, combining the stipulated propagator properties, one starts by equivalently writing where the last identity follows since using the identity between the retarded (+) and advanced (-) Green functions [64] ( combined with the propagator-Green function relationship (130), while supplemented with the advanced propagator version Now, the propagators from (153) may be written immediately under the general form (152) contributing in rewriting (153) as (160) With expression (160) the propagator (152) is fully expressed in terms of classical action as or in the more appealing form usually referred to as the Van Vleck-Pauli-Morette formula, emphasizing on the importance of solving the classical problem for a given canonical Lagrangean [60,66].
However, the path integral solution (161b) has to be used with two amendments: o the procedure is valid only when the quantity (158), here rewritten in the spirit of (161b) as , performed respecting one end-point coordinate remains linear in the other space (end-point) coordinate b x , so that the identity (159) holds; this is true for the quadratic Lagrangeans of type (146) but not when higher orders are involved, when the previously stipulated Fourier analysis has to be undertaken (one such case will be in foregoing sections presented). o In the case the formula (161b) is applicable, i.e., when previous condition are fulfilled, the obtained result has to be still verified in recovering the delta-Dirac function by the limit in accordance with the implemented recipe, see Equation (153); usually this step is providing additional phase correction to the solution (161b). The present algorithm is in next exemplified on two paradigmatic quantum problems: the free motion and the motion under harmonic oscillator influence. In each case the knowledge of the classical action will almost solve the entire path integral problem.

Path Integral for Free Particle
Given a free particle with the Lagrangean it leads by means of Euler-Lagrange equation to the classical (Newtonian) motion fulfilling the boundary conditions a a cl being these endpoints the states where the system is observable, i.e., where the quantum fluctuations vanishes, see Equation (145c) and Figure 2.
Replacing solution (165b) back in Lagrangean (163) the classical action is immediately found Next, the quantity (158) is firstly evaluated in the spirit of (161b) as and recognized as linear in the other end-point space coordinate b x . Thus, the formula (161b) may be applied, with the actual yield Finally, the result (168) has to be arranged so that to satisfy the limit (162) as well. For that we use the delta-Dirac representation Comparison between (168) and (169) leads with identification thus correcting the factor of (168) towards the correct limiting path integral solution Remarkably, this solution is indeed identical with the Green function of the free particle, up to the complex factor of (130), thus confirming the reliability of the path integral approach. Moreover, beside of its foreground character in quantum mechanics, the present path integral of the free particle can be further used in regaining the energy quantification of free electrons in solid state (motion within the infinite high box) as well as the Bohr quantification for the continuous deformation of the path on the circle [39,65].
Yet, these cases appeal the spectral representation of the quantum propagators and will not be treated here, being more suited for a dedicated monograph [67].

Path Integral for Harmonic Oscillator
The characteristic Lagrangean of the harmonic oscillator ( ) provides, when considered in the Euler-Lagrange equation (164), the classical equation of motion with the well known solution specialized for the end-point events of motion as ( ) In the same way as done for the free motion, see solution (165b), worth rewritten the actual classical solution (173b) in terms of relations (174), for instance as or similarly as: On the other hand the classical action of the Lagrangean (172) looks like Now, in order to have classical action in terms of only space-time coordinate of the ending points, one has to replace the end-point velocities in (176) by the aid of relations (175a) and (175b) in which the current time is taken as the b t t = and a t t = , respectively; thus we firstly get then we form the working products to finally replace them in expression (176) to obtain the computed classical action Note that the correctness of Equation (179) may also be checked by imposing the limit 0 → ω in which case the previous free motion has to be recovered; indeed by employing the consecrated limit Such a kind of check is most useful and has to hold also for the quantum propagator as a whole. Going to determine it one has to reconsider the classical action (179) so that the quantity (158) is directly evaluated in the spirit of (161b) as thus again encountering the linearity case in the other end-point coordinate b x . Being this the fortunate situation in which the previous expression (161b) for path integral computation may be applied with the harmonic oscillator result Yet, as above was the case for the classical action itself, also the pre-exponential quantum fluctuation factor of (183a) has to overlap with that appearing in the path integral of free motion of Equation (171) for the limit 0 → ω ( Thus we have to adjust the propagator (183) with the exponential pre-factor corrected with the complex factor " i " This is the sought propagator of the (electronic) motion under the harmonic oscillating potential, computed by means of path integral; it provides the canonical density to be implemented in the DFT algorithm (128), (129) Yet, for practical implementations, the passage from quantum mechanics (QM) to quantum statistics (QS) is to be considered based on the Wick transformation (10) here rewritten as providing the Euler trigonometric to hyperbolic function conversions (by analytical continuations) finally displaying for density (186) the counterpart formulation The uni-particle (electronic) density (189) is then used for computing the harmonic oscillator partition function Now, using the "double angle" formula 1 sinh 2 1 cosh 2 sinh cosh 2 cosh 2 2 2 2 Remarkably, the result (192) recovers also the energy quantification of the quantum motion under the harmonic oscillator influence, as seen by the successive transformations there follows immediately the harmonic oscillator energy quantification 194) in perfect agreement with the consecrated expression.
The results of these two sections suggest the following rules for using path integrals propagator for density computations: o The reliable application of the density computation upon the partition function algorithm, see Equations (128) and (129) with the complex factor " i " included, as confirmed by both the free and harmonic oscillator quantum motions; it may be used for linear classical actions in one of the end-point space coordinates upon derivation respecting the other one; yet the formula (195) should be always checked for fulfilling the limiting (162) delta-Dirac function for simultaneous events for any applied potential. Nevertheless, recognizing the major role the classical action plays in the path integral representation of the quantum propagation (and propagator), the question whether it is possible to consider the semi-classical expansion of the propagator in general case, without being under any constraint except the semiclassical (higher temperatures) limit itself, naturally arises. Such an approach is exposed and its reliability tested in the next sections.

Semiclassical Expansion
Semiclassical derivation of the evolution amplitude employs some of the previously Feynman path integral ideas refined due to the works of Kleinert and collaborators [39,41,[68][69][70][71]. They are bellow summarized.
o The real time dependency is "rotated" into the imaginary time according with the Wick transformation (10). o The quantum paths of (145a) are re-parameterized as where the classical path of (145a) is replaced by the fixed (non time-dependent) average while the fluctuation path ( ) τ η remains to carry the whole path integral information, yet being departed at the end of integration frontier from previously Dirichlet boundary conditions (145c), where it vanished at the domain frontiers, to the actual different endpoint values In these conditions the quantum statistical path integral representation of quantum propagator becomes since we immediately noted the immediate transformations based on the above (197)-(200) parameterization. It should be pointed out that the used re-parameterization is not modifying the value of the path integral but is intended to better visualize its properties, towards evaluating it. As such, from expression (201) it now appears clearer than before that for the systems governed by smooth potentials, the series expansion may be applied respecting the path fluctuation, here in the second order truncation where the covariant notation for products was assumed for maintaining the generality of the Ddimensioned approach. This way, a (truncated) series of path integral evolution amplitude of Equation as being driven by the quantum fluctuation' various orders contributions, up to the second order. This is a natural approach since the very quantum nature of the path integral is given by the quantum fluctuations themselves, from where the systematic approximations of path integrals over the quantum fluctuations. The series is known as the semiclassical expansion since is formally done in the "powers of h -Planck". Now, looking on Equation (204) where we have used (196b) along identifying the semiclassical factor . However, the expression (205a) may be further formally cast by introducing the so called free imaginary time amplitude readily given by the free-propagator solution (171) accommodated by the present statistical and boundary transformations Equations (196b) and (199) to the form while having the normalization role for averaging the semiclassical factor contribution Therefore, the semiclassical form of path integral representation of evolution amplitude looks like The remaining problem is that of expressing the averaged values of the fluctuation paths in single or multiple time connection, i.e., From the heuristic point of view it is normal to arrive at the form (208) because it tells us that the quantum fluctuation is firstly averaged along the quantum evolution and then averaged by time in order that the evolution amplitude is determined. Observe also that the present semiclassical approach is not using the previously employed properties of the classical action, avoiding therefore the limitation of the derivative behavior at edge of the space domain of integration, while posing now the limitation in what respect the quantum fluctuation power. It is also useful to remark that the present semiclassical approach may use the interplay between the previous solved free-and-harmonic quantum motions, since the path integral (206a) may equally be regarded as the free motion of the quantum fluctuations (naturally since they are not known a priori or with some possibility of instantaneously observation); at the same time, if one formally counts the kinetic term as the perturbative (a.k.a. fluctuation) oscillatory motion a more complex picture of quantum fluctuation is obtained; in conclusion, quantum fluctuating paths may be (or should be) treated as being a kind of harmonically free motion: harmonic since as fluctuations may be expanded in Fourier series (as originally perceived by Feynman), but also free since their unknown of instantaneous feature. Therefore, an appropriate use of both these manifestations will conduct to the reliable path integral representation. Yet, since we have already used the free-motion character of fluctuation paths, the harmonic one is next entering the analysis.

Connected Correlation Functions
For calculating the average of quantum fluctuation paths one has to understand their inner nature: in order reconciliation of free and harmonic features be achieved the so called quantum current ( ) τ j is introduced (and presumed to appear in reality too as causing/driving the quantum fluctuations), so that the propagator of this current, known as the generating functional, is formed [39,66] ( ) with which help one can recognized the equivalence in accordance with general definition (207). One can nevertheless see that the quantum current appearance in (211) is under the perturbation form, so that it readily accounts for the deviation from the free fluctuation motion towards the harmonically one. Therefore, although general correlation definition may be advanced by the ordering rule the problem of practically evaluation still remains. Aiming for solving it one observes that the form (212) is analogous with the partition function based electronic density, see Equation (12) for instance; consequently, the alternative formulation looks at the canonical (N = 1, mono-particle) level like where, now, the quantity [ ] with + S being the Euclidian action, see Equation (121c), while space and time slicing intervals are those introduced by (81) and (90), respectively. The disconnected character of correlations (213)-(214) may be overcome remembering that the logarithm of the partition function provides the thermodynamic free energy, see relation (67) Nevertheless, aiming to have a better "feeling" on how the connected and disconnected correlation (fluctuation) functions (217) and (213) are linked, let's start evaluating some orders of them.
As such, absorbing the constants in the involved functionals, the first order of (213) reads correlation of (213) we successively have as the single connected path (remember that the fluctuation was already averaged out) that is just the classical path connecting the ending points of the quantum evolution, see Figure 2. Now, going to the second order of correlation of (213) one has a result that can be wisely rearranged as or, even more practically for our purpose, as In similar manner, while applying a kind of recursive rule, sometimes denoted as the cluster decomposition or cumulant expansion while involving the pair-wise (Wick) decomposition of the n-points correlated function one can easily obtain the higher orders of correlations, however observing that all connected orders of events reduce to the combinations of pair-connected events. For instance, we get for the third order fluctuations the average contribution or with more terms involved for the fourth order τ τ τ τ τ τ τ τ τ τ τ τ con con con con con con Next, having these examples in hand, one tries to re-deriving them by an appropriate generating functional (210) worked with the connected function definition (212). At this moment one uses the previously emphasized "free-harmonic motion" dual nature of fluctuation paths -see Equation (209), to reconsider the free imaginary time amplitude (206a) contribution with harmonic Euclidian action of fluctuations while fulfilling (at the end of calculation) the constraint We like to rearrange the action (222b) so that the quantum current contribution to clearly appear; in achieving this one firstly rewrites it by performing the integration by parts where we have recognized the harmonic differential operator The form (223) is very useful through employing the Green equation of harmonic motion for the integral property to perform the path shifting of fluctuations by the transformation so that the prescribed action of (223) takes the form under the assumption the physical integration interval ( ) b a τ τ , assimilates the entirely evolution universe of the concerned problem, i.e., the mathematical interval ( ) , so that the delta-Dirac integration property (225b) is consistent. Also note that in expression (228) since the statistical Green function should come from its associate original real time quantum mechanically problem, see below, it tracks also the temporal Wick "rotation i t / τ = " in the integration measure, explaining therefore the complex factors in the last term of (228).
With these, the harmonic fluctuation action of (228) may be reconsidered with the working form such that it can be further rearranged so that the free terms action to appear distinctively under the with the free Euclidian action recovered though considering on expression (229) the reverse integration by parts -as unfolded from (222b) to (223). Finally, back with the identifications in the action (230), there follows that its last two terms are no longer displaying quantum fluctuations upon integration since they were comprised under averaged forms, see Equations (218) and (219c), respectively; thus they release for the searched current-dependent amplitude of (222c) the actual solution which ultimately simply re-writes like the current dependent propagator amplitude as a wave-perturbative form of the free fluctuation amplitude (206a), being intermediated by the harmonic towards free limiting motion of quantum fluctuations. Let's further comment that the actual form (232b) generalizes the previously "guessed" form (210), which provided the first order fluctuation correlation, however having in addition the power to recover all other superior orders of correlation, for instance those given by Equations (219c) and (221), by successive application of the formula (212); this way, the transformation factor m / h in (231b) is as well justified. With these the connected correlation function algorithm was proofed in detail, being at disposition to be implemented for whatever order of semiclassical expansion of the path integral evolution amplitude (208); as exemplification, the next section will expose the analytic solution for the second order case.

Classical Fluctuation Path and Connected Green Function
We have already seen that aiming to evaluate any of the above connected correlation functions one imperatively needs to know the analytical forms of classical fluctuation path and for the connected Green function -identified from (231b) as while depending on its turn of the knowledge of the Green function for the harmonic oscillator problem, Equation (225a) with (224), to be finally specialized towards the "free harmonic" limit 0 → ω . Therefore, with the quantities of Equations (233) any semiclassical problem can be solved analytically. Yet, for the quantum objects in question the computing procedure consists by three major stages, as unfolded in the sequel.
(i) Solving the associate real time harmonic problem; (ii) Rotating the solution into imaginary time picture; (iii) Taking the "free harmonic limit" 0 → ω .

Calculation of Classical Fluctuation Path
As discussed above, the classical path for quantum fluctuation will not be written directly from the ordinary path free motion (Section 3.3.2) but using the similar result of harmonic motion (Section 3.3.3) upon which the free-harmonic condition 0 → ω will be imposed; actually, the procedure is unfolded as follows.
The result (175a) is combined with (177a) to provide the real time classical path The real to imaginary time rotation is performed on the result (234) according with the Wick rule prescription of (196a), being this equivalently of directly rewriting of expression (234) replacing the trigonometric functions by their hyperbolic counterparts, according with the previously explained conversion, see Equation (188a) The "free-harmonic" ( 0 → ω ) limit is performed upon the expression (235) through employing the ordinary hyperbolic limit which evidently does the same job as the classical free-motion result of (165b), although not identical, since derived from a generalized perspective here.
The result (237) is implemented in the formula (197) to finally produce the classical fluctuation path which takes even the simpler form when rewritten within the thermodynamic picture

Calculation of the Connected Green Function
Now, going to the evaluation of the expression (233b) we need the Green function of the harmonic oscillator from the equation (225a); when rewritten in real time picture it has the advantage of having the frontier values fixed by the Dirichlet boundary conditions ( ) in the same manner as the fluctuation paths in real time are set to vanish at the endpoint frontier, see Equation (145c). Such double boundary condition fixes the type of solution as being of the double trigonometric form In the same manner the temporal alternative ordering problem of (240) with the same constant as for the solution (240c) since recognizing that both formally belong to the same homogeneous equation of type Now, looking for appropriate identification in the inhomogeneous equation one notes that its left side is formed from the difference of the first derivatives of the solutions (240c) and (241c) approaching each other for the concerned times which, through comparing with the right side first term of (243) gives the searched constant It leaves with the real time Green function solution of the harmonic oscillator that nevertheless combines both above solution with the help of Heaviside step-function Next, as previously done with the fluctuation paths, the change to the imaginary time picture is done automatically through trigonometric-to-hyperbolic recipe (188a) to give noting that in the course of transformation the factor was tacitly absorbed, with the parenthesis complex indices coming from the trigonometric to hyperbolic rotation (188a), while the outside index complex assures the equivalence of Green function contribution for the canonic-to-Euclidian path integrals action exponents' transformation Expression (248) is finally employed to the "free harmonic" limit (236) providing the result which being free of harmonic influence it remains identically also from the quantity ( ) it has to be converted into the searched connected Green function (233b), leaving with the time imaginary form or with its equivalent statistical one when the thermodynamical picture (239) is considered.

Second Order Semiclassical Propagator, Partition Function and Density
Returning to evaluate the second order truncated expansion (208) one needs the evaluation of the quantities and of their integration. Given the previous discussions, see for instance the Equations (218) and (219c), one immediately has With the help of expression (238) the first order averaged fluctuation integral appearing on (208) becomes obtained by replacing into the expression (252b) the classical fluctuation paths and the connected Green function components, Equations (238a) and (251a), respectively. Now, the second order averaged fluctuation integrals are computed as following: • At coincident times or re-written as ( ) ( ) within the thermodynamical environment given by Equation (239).
• At different times with its quantum thermodynamical counterpart Finally, while replacing the values of Equations (253), (255b), and (256b) back in the second order truncated semiclassical expression of imaginary time amplitude (208) one gets the analytical result Note that the expression (257) plays the role of the semiclassical canonical density in PI-DFT algorithm given by Equations (128) and (129) At this point, the expression (260) may be elegantly transformed through considering the Gauss theorem of integrated divergence that written in a general D-dimensional case leaves with the useful differential relationship helping in rearranging the partition function (260) firstly as and finally, after the exponential resume, equivalently as In the same manner also the higher orders of semiclassical expansion of density matrix (204) or (208) can be constructed by following the cumulant expansion (220), its fluctuation path and connected Green function components, as given by Equations (238) and (251), respectively, towards producing the analytical canonical density, the partition function and finally the many-body density to be used in density functional theory and of its (chemical) applications [71]. Such an application is to be in next presented for electronegativity and chemical hardness indices' computations.

Fourth Order Semiclassical Electronegativity and Chemical Hardness
Here, we assess electronegativity (EN) as the convolution of the imaginary time conditional probability ( ) , with the radial valence shell potential ) (r V [72] ( )dr so representing the power of the entire atom (nucleus + core electrons + valence shell) to attract electrons of the outer shell (fixed by radius r ) to its center ( 0 = r ). This way, the current EN definition is refined to account for the whole stability of the atom with its electronic and nuclear subsystems.
Having an analytical EN quantum formulation, the chemical hardness, η, its natural companion, is re-expressed from Equation (13b) under the explicitly working form [73] playing a major role in establishing the main chemical principles of reactivity towards stability: the hard-and-soft-acids-and-bases (HSAB) and the maximum hardness (MH) [74][75][76][77][78].
The general radial one-dimensional probability amplitude connecting the space-time events (r a , 0) and (r b , τ) for electronic evolution amplitude in an atom can be derived from the semiclassical expansion, by extending the expansion (208) up to the fourth order through including the terms of type (221), whose calculation leads to the result [71] ( ) where eff Z stands for the Slater effective atomic number specific for the multi-electronic atoms, being derived from the standard atomic number Z by subtracting the shielding effects of the inner electrons [79]. Nevertheless, worth mentioning that the usual negative sign in attractive potentials was formally abolished because: we retain the positive values of electronegativity (263) since EN is evaluated as a stability measure of such nuclear-electronic system; the sign is in accordance with the electric field orientation that drives the sense of the electronic conditional probability of the imaginary evolution amplitude evaluated from the center of atom ( 0 = a r ) to the current valence shell radius ( b r ).
Therefore, the electronegativity can be seen also as power of holding electrons in the valence shell reciprocal to that exercised upon them from the center of atom. This way, the present EN definition and equation stand for the reconciliation of the two opposite phenomena acting upon the valence electrons: attraction to nucleus and repulsion from the other atomic inner electrons.
Next, within the Bohr description of the electrons moving in a central potential [80], while adopting the atomic units, However, for keeping the analyticity of the present approach, the computation of the integral (263) with the replacements (265)-(268) may use the saddle-point recipe, very well accommodated for the present semiclassical context; thus we implement the approximation rule [81] corresponding to the valence shell saddle radius expressed out by the optimization condition According with the exposed strategy, one finds the fourth order semiclassical expansions for electronegativity is given [   The numerical fourth order semiclassical electronegativity and chemical hardness atomic scales are reported in Table 2 and represented in Figure 3, where the comparison with the finite-difference counterparts of Table 1 was also emphasized. The striking difference in terms of orders of magnitudes observed between elements down groups is the main characteristic of the actual atomic scales of electronegativity and chemical hardness; however, due to the fact the actual definition of electronegativity and chemical hardness reflects the holding power with which the whole atom attracts valence electrons to its center -this is not a surprising behavior.
It is therefore natural to observe that as the atom is richer in core electrons down groups lesser is the attractive force on the outer electrons from the center of the atom. In this regard, the actual scales mirror the atomic stability of the valence shell at the best. Nevertheless, a better regularization of their increasing trend along periods it is observed in Figure 3 for both semiclassical electronegativity and chemical hardness fourth order scales, a feature more apparent for the actual chemical hardness scale. Moreover, since the chemical hardness controls the secondary order effects throughout its definition as the derivative of the electronegativity a phenomenological rule would demand to have lower values than that of the associated electronegativity.  Table 2 respecting their finite-difference (FD) counterparts of  [83]: • the atoms N, O, F, Ne, and He have the highest electronegativities among the main groups; • the electronegativity of N is by far greater than that of Cl -a situation that is not met in the finite-difference approach; • the Silicium rule demanding that most metals to have EN values less than or equal to that of Si, is as well widely satisfied; • the metalloid band (B, Si, Ge, As, Sb, Te) clearly separates the metals by nonmetals' EN values; • along periods the highest EN values belong to the noble elements -a rule not fulfilled by the couples (Cl, Ar), (Br, Kr), and (I, Xe) within the finite difference representation, see Table 1; • the recorded electronegativity values of the chalcogens (O, S, Se, Te) reveal great distinction between the chemistry of oxygen and the rest elements of VIA group; • the transitional metals are grouped in a distinct contracted region of EN values -this way closely emphasizing on the d-orbitals effects, a criteria almost not fulfilled by the finitedifference scheme, see Table 1 and the Figure 3. Finally, we have to point out that the systematic decrease of orders of magnitude of electronegativity and hardness semiclassical scales of Table 2 and Figure 3 has a fundamental consequence, namely it stands as the computational proof that the electronegativity and hardness behave like pure quantum/structural indices. As such, they are not manifesting with the same intensity among all elements of the Periodic Table by having values that tend to considerably diminish as the frontier electrons are farer and feel less and less the quantum influence (potential and force) of the nucleus and of the core electrons. These results are in accordance with the electronic localization principles in an atom [3,84].

Effective Classical Partition Function
As previously shown, see Equations (121c) and (208), for instance, considering the path integral propagator that underlies the canonical density in the quantum statistical algorithm, see Equations (87)-(89), accounts for the quantum effects (fluctuations) induced on single particle paths by the presence of an external potential, while being analytically computed by averaging these over all possible configurations. Yet, one could observe that for periodic paths, i.e., when the final and initial space-points coincide b a x x = (274) the particle travels in very short time not far away from the initial position and then is back on the initial point; such picture has the physical measurable consequence a particle is observed on the initial point, i.e., it is found on a stationary state/orbit, while the quantum fluctuations are oscillating around the equilibrium (initial = final) space-point. Even clearer, the situation corresponds to the classical picture in which a particle behaves, being accommodated in an equilibrium state/stationary orbit under external potential influence. This means that the external influence itself is observable in (initial = final) concerned/measured state, thus being no longer a path parameterized function, but a constant: Therefore, the associated periodic propagator (density matrix) becomes where the recognized path integral of free motion was solved by plugging into its quantum mechanical solution (171) the present conditions (274) and (196b), for accounting of the path periodicity and quantum statistics, respectively. At the same time there is clear that the periodic path condition (274) is not arbitrarily but a compulsory step since characteristic in passing from density matrix to partition function and then to the real (measurable or workable) canonical and N-particle density, according with the density matrix algorithm (87)- (89). Therefore, the resulting partition function built from the un-normalized canonical density (276) assumes the simple form is approximated with unity in expression (276), thus with an error of 40% at the maximum displacement of (278) value; as the classical displacement (278) tends to zero as the expressions (276) and (277)  With these considerations there appears as natural the generalization of the classical partition form (277) into the more comprehensive one known as the effective classical partition function [85][86][87][88][89][90] [ ] with the integration variable defined as the thermal average of the periodic quantum paths sometimes called as the Feynman centroid, while the notation is to be right bellow justified. Moreover, the search for the best approximation of effective-classical partition function (280) will be conducted as such the quantum fluctuations be not dependent on the classical displacement (278), abstracted from the free motion, but being driven by the quantum harmonic oscillations -through they constitute a generalization of the free motion itself, see for instance the equivalence of classical paths or propagators of free with harmonic motion in the zero-frequency limit, see the Section 3.3.3.
However, the periodicity condition (274) for paths is to be maintained and properly implemented in approximating the effective-classical partition function (281) being, nevertheless, closely and powerfully related with the quantum beloved concept of stationary orbits defined/described by periodic quantum waves/paths. This way, the effective-classical path integral approach appears as the true quantum justification of the quantum atom and of the quantum stabilization of matter in general, providing reliable results without involving observables or operators relaying on special quantum postulates other than the variational principles -with universal (classical or quantum) value.

Matsubara Frequencies and the Quantum Periodic Paths
As always done when a new type of path integral is under consideration the reconsideration of the quantum paths, and in fact the quantum fluctuations, is undertaken so that facilitating the best way for solving it. Yet, this time due to the periodicity condition of paths the propagator is hidden by the associated partition function. Therefore, the optimum approximation for the effective classical potential in (281) will provide the periodic evolution amplitude as well, i.e., the un-normalized density, which by normalization with partition function will lead with the searched canonical density counterpart.
Going to characterize the periodic paths, they will be seen as the Fourier series [34,39,41,91,92] ( ) ( ) in terms of the so called Matsubara frequencies m ω ; they are explicitly found through specializing the condition (274) into the actual statistical one, see Figure 4 b a resulting in the equality   with the 0 th terms viewed more than the "zero-oscillating" or free motion path but the thermal averaged path over entire quantum paths (282a) However, beside revealing the integration variable of the classical partition function (281a) as being of averaged nature, the result (282e) emphasizes on the actual periodic path decomposition (282d) integral featuring another level for parameterization of quantum paths that goes beyond characterizing them as quantum fluctuations around classical motion; they are here constructed as periodic oscillations (back and forth -due to their complex form, in analogy with conjugated plane waves traveling in opposite directions) around the averaged path value (interpreted as thermic average, or, more plastic, as centroid of the quantum fluctuations themselves). Therefore, with such path parameterization perspective the present level seems involving quite complex quantum phenomenology to be further enriched in the sections to follow.

Matsubara Harmonic Partition Function
The quantum path decomposition (282d) imposes the factorization of the Feynman path integral measure (120) accordingly with the integration constants 0 C and m C to be determined from identifying the known partition function of the harmonic oscillator, Equation (192), with the path integral representation of the counterpart partition function written by the measure (287) and paths (282) To this end, let's begin with the computation of the kinetic term appearing under the integral (288a) while for the harmonic contribution we similarly get  Going back, once calculated, the Riemann series (296) is replaced in (295) which, at its turn, is employed for the integral evaluation for ending with unfolding of the searched function (293). These steps will be systematically exposed in next sections.

The Generalized Riemann' Series
Computation of the generalized Riemann series (296) in terms of the extended series ( ) o Computing the integral under the sum of (299) by the complex integration, according with the contours of integration identified in Figure 5 around the poles α i q ± = throughout applying the residues' theorem   o Insertion of these integration results in the expression (299) is done by attributing to each contour and integration the (series) summing range according with the constraints of (301a) to successively yield We have now clarified all prerequisites to readily compute the Matsubara harmonic partition function (292), used as a tool to find out the Matsubara normalization of periodic path integrals. This will be addressed in the sequel.

Periodic Path Integral Measure
The Matsubara harmonic partition function algorithm (292)-(297) may be now unfolded successively as: o Computing the function (295) by inserting the above Riemann generalized series (304): o Evaluating the function (294) by the aid of (297) rule through considering the variable change Note that the measure given in (316) is rather universal for periodic paths, while the involvement of the harmonic oscillator was only a tool (and always an inspiring exercise) for determining it since the complete quantum and statistical solution at hand. As such, once more, the harmonic motion proofs its versatile properties respecting the fluctuation over -or perturbation of -the free motion by modeling the quantum displacements from classical equilibrium or observed path.

Feynman-Kleinert Partition Function
Being equipped with the periodic path integral technique we can present one of the most efficient ways for approximate the effective-classical partition function (281a); it starts with the general path integral Since the analytical solution for expression (317) is hard to be conceived for an unspecified potential form, it may be eventually reformulated in a workable from by involving another partition function, the so called Feynman-Kleinert partition function FK Z and its special average recipe In Equations (318a) and (318b) the Feynman-Kleinert partition function writes as [85] with the working action ansatz assures the global optimization for the action, and implicitly for the Feynman-Kleinert partition function, so approaching at the best the exact partition function (317) and its associate total ground state energy of the system given by the free energy In fact, the Feynman-Kleinert action (319b) is to be involved in two-fold optimization algorithm for providing the best approximation of the partition function (317). This will favor a close analogy with the double search for electronic density, in density functional theory (DFT), as will be discussed later.
Yet, the Feynman-Kleinert partition function is to be unfolded within the actual periodic path integral representation ( ) ( ) ( ) to be optimized respecting its harmonic frequency (for equilibrium optimization) and for ground state perturbation (optimization) in what follows.

Feynman-Kleinert Optimum Potential
The optimization of the Feynman-Kleinert partition function (322) is performed employing the Jensen-Peierls inequality whose the phenomenological proof is given in the Figure 6, on the partition function relationship (318a) leading to the lower bounded partition function or, by calling the Equation (320), to the higher bounded free energy When rewritten the last inequality with the help of Euclidian actions for general partition function and the Feynman-Kleinert specialization, Equations (317) and (319), respectively, one notes the cancellation of the kinetic (free motion) terms, while the resulting expression provides the searched variational architecture where all involved terms combine the external, perturbation and quantum fluctuation influences.
Having the variational problem formulated it remains to individually compute the terms appearing in the Feynman-Kleinert average (328), by using the associate definition (318b) with the action (319b).
Going to evaluate the most general term containing the external potential average, we have in the first instance its periodic path integral representation where the Fourier k-(wave vector) representation was implemented for external potential so that the quantum path to explicitly appear in evaluation; this technique had helped for performing the quadratic completion of paths in the view of harmonic-like integration of type (290)  which can be immediately recognized as directly related with generalized Riemann series (304), thus having the form [85] ( ) ( ) Expression (329) can be even more simplified when solving out the k-integral by back considering the Fourier transformation for the potential and then proceeding with the quadratic completion toward the Poisson standard integration The potential (331) is known as the smeared out potential and has a major role in explaining the quantum stabilization of matter, as will be largely discussed in the next section. For the moment it is regarded jus as the integral transformation of the original applied potential by convolution with a Gaussian packet with the width ) ( 0 2 x a that accounts for the existing quantum fluctuation in the system. Nevertheless, Equation (331) leaves the Feynman-Kleinert average of external potential with the result [85] For the rest of the averaged terms in (328) the evaluations are considerably easier since for each of them we have only to compute their smeared out version (331), with the yield for the trial harmonic and respectively for the Feynman-Kleinert trial function from where the first stage of variational algorithm is fulfilled by the obvious choice With Equation (336) the Feynman-Kleinert potential (323) now displays as [85] ( ) ( ) ( ) ( ) There remains only to finally optimize the explicit potential (337) for the harmonic (trial) frequency assuring therefore the equilibrium of the gained lowest approximation of the ground state for the concerned system. This is simply achieved through the chain derivative Nevertheless, through observing the huge role both the smeared out potential (331) and the fluctuation width (330) play in deriving the approximated equilibrium ground state they deserve be further analyzed and commented in relation with matter stability.

Quantum Smeared Effects and the Stability of Matter
The intriguing role the smeared potential in special and the smearing effect in general play in optimization of the total energy and partition function of a quantum system opens the possibility analyzing the "smearing" phenomenon of the quantum fluctuation in a more fundamental way.
I. Firstly, it was noted that the smearing potential (331) appears as a Gaussian convolution of the applied potential, although modeling the evolution of a wave-packet under that potential influence; in other terms, it appears the fundamental question whether the Gaussian and wave function "kernels" behave in similar way throughout the smearing effect of quantum fluctuations; analytically, one likes to see whether the next smearing average equality readily holds In order to check (342) one separately computes each of its sides by the aid of k-form of (331) and successively gets the smearing average for the wave-function and for the Gaussian packet Now, for closely comparison of the expressions (343a) and (343b) the most elegant way is to make once more recourse to the smearing procedure, this time referring both to the entire paths and Feynman centroid; to this end, the previous result (333) is here used in the variant: ( ) Since the difference between these expressions is numerically proportionally with the factor ( ) 029 they can be considered as identical in quantum smearing effects and Equation (342) as valid. Yet, the quantum identity between the plane-wave and Gaussian packet has profound quantum implication, while revealing for instance the de Broglie -Born identity in Gaussian normalization of the de Broglie moving wave-packet. It may express as well the observational Gaussian character of the wave-function evolution in Hilbert space. Finally, and very important, it leads with explanation of the Bohr first postulate, i.e., it is able to explain the stationary wave on orbits under singular (Coulombic) potential thus explaining the matter stabilization on rigorous quantum base, rather than to admit it by the power of a postulate. This is to be proved next [39,85].
II. Let's consider a quantum system evolving under the influence of the Yukawa potential, as a generalization of the Coulomb interaction, available also for the sub-nuclear world , which goes to the celebrated Hydrogen Coulomb central potential in the limit ( ) thus assuring (and explaining) why the atomic electron(s) do not fall onto the nucleus. Therefore, the smearing procedure plays a kind of renormalization role in transforming singular potential in finite interactions by means of quantum fluctuation effects. Such picture strongly advocates for the powerful path integral formalism in general and of that of Feynman-Kleinert in special for explicitly accounting of the fluctuation width in optimizing the quantum equilibrium states. Nevertheless, worth particularizing the Feynman-Kleinert formalism to the ground and excited states for better capture its reliability and limits.

Ground State (β→∞, T→0K) Case
The basic ground state conditions in terms of thermodynamic factor ( β ) or the temperature (T ): aims to bring the Feynman-Kleinert formalism, through its working potential (337), closer to the ground state as usually provided by the consecrated quantum variational principle; for this purpose it will be firstly specialized within the general limit (357) and then tested for the paradigmatic Hydrogen ground state case for investigating upon the accuracy of the formalism itself. As such, the first component of the Feynman-Kleinert potential (337) has the ground state limit recovering the ground state of harmonic motion for the trial fluctuations, while the ground state of the fluctuation width (330b) reads as from where also the asymptotic trial fluctuations' frequency springs as Next, by combining relations (358) with the working general effective-classical approximation potential (337) the general ground state limit looks like The ground state smeared out potential remains to be individuated in solving the associated ground state problem.
Very interestingly, the expression (359) entirely corresponds to the smeared out effect applied on the ordinary quantum Hamiltonian ( ) as one can immediately check out through applying the general smearing averaging definition (331) on it The identity between expressions (359) and (361) leaves with the important confirmation that the smearing operation produces in fact the average of quantum fluctuation for the ground state equilibrium. For the Coulomb interaction, say for the Hydrogen atom, either of above expressions produces the working form where the 3D version of the kinetic term of (361) was here considered aside the smearing out potential in the origin (356b); they turn out the ready form for ordinary minimization respecting the fluctuation width thus producing only a 6% error in predicting the radial localization in stabilizing the electronic ground state orbit, i.e., being closer to the nucleus respecting the exact Bohr-Schrödinger solution. However, for predicted approximated ground state energy the error is a bit higher due to the specific dependency predicting the spectral localization with an error about 16% higher than the exact ground state of Hydrogen atom.
Nevertheless, besides the approximated character of the formalism, the Feynman-Kleinert approach adapts very well to the singular potential, having the advantage of being compatible with a wide class of electronic potentials in atoms and molecules; moreover, it can be particularized to the ground state in the same degree it accounts for higher temperature cases, the other extreme of thermodynamic limit -see the next section, spanning this way in principle the entire range of statistical systems at equilibrium; note that such "universal" quantum statistical picture of equilibrium is hard to found in the quantum theory, at the same level of elegance, analyticity and complex ideas [93][94][95][96][97].

Excited State (β→0, T→∞) Case. Wigner Expansion
As before, the components of the Feynman-Kleinert potential (337) are to be now evaluated in the excited states or the valence limit This limit is applied on terms of Feynman-Kleinert potential reduced to the hyperbolic functions on which the approximations of type (306) may be applied.
With this recipe we firstly evaluate where, in the last expression, the hyperbolic cosecant was approximated, in the spirit of Equation (306), with the form With the partial limits (371) the expansion (370) becomes providing on its turn the harmonic term approximation ( ) ( ) with the same first β -order dependency as the limit (368) of the fluctuation width.
The remaining term for evaluation in higher temperature limit is the smeared potential (331); for it the next change of variable will be done by considering allowing for the successively series formulation known as the Wigner expansion [98] ( ) The expression (376) was obtained due to the small fluctuation width at higher temperatures, see the limit (368) and Feynman's discussion in Section 5.1 when introducing the idea of effectiveclassical potential and partition function; it features the small perturbation in terms of fluctuation width around the applied potential "centered" on the Feynman centroid, therefore behaving as a sort of semiclassical expansion.
Indeed, by recalling the trial fluctuation frequency optimal definition (341) it specializes in the high temperature limit potential (376) to the working expression Finally, by plugging the optimum frequency (377) into the limit (374), and along the limits (368) and (376), back in Feynman-Kleinert potential (337) it acquires the high temperature or the valence form Most remarkably, the last form is nothing than the semiclassical potential appearing in the second order partition function (262b), thus providing the identical Feynman-Kleinert partition function (322) From such identity one may re-affirm the important conceptual achievement according which the Feynman centroid (281b) for periodic path approach corresponds with the end-point coordinate average (198) in semiclassical expansion.
Moreover, there is clear that the higher temperature or the valence limit of the Feynman-Kleinert periodic path integral approach regains the semiclassical expansion of the non periodic paths of a quantum particle; they become nevertheless periodic at higher temperatures due to higher oscillations around the equilibrium, while possessing not sufficient kinetic energy to break the equilibrium by traveling too far away.
We are now fully convinced that Feynman-Kleinert path integral formulation works fine either at low and higher temperature limits, while recovering both the (Hydrogen) ground state and the semiclassical valence expansion with reliable fidelity, respectively. Nevertheless, it remains to stress on further connection with the electronic density and consequently with the density functional theory towards the quantum chemical applications. These aspects will be addressed in the next section. With these, the working form for the canonical electronic density can be written in terms of the Feynman centroid (281b) by the expression

Path Integral Connection with Density Functional Theory
satisfying the centroid normalization ( ) or, more sophisticatedly, in terms of the quantum path (282a) with the respective normalization to unity Remarkably, the Feynman-Kleinert variational algorithm in path integrals may be viewed as providing the calculation of the effective electronic density by constructing the constraint-searched partition function picture as the Levy constraint-search formalism [99] does in seeking the electronic density from the trial wave function. More clearly, Levy's recipe prescribes that the ground state energy minimization scheme (within the second Hohenberg-Kohn theorem [1]) involves, in fact, two steps: one over all wave functions ( Ψ ) that give the same density (the inner minimization step) followed by the minimization throughout all density classes (the outer minimization step) through recalling the density functional basic functionals (3)- (5).
Such equivalence between the path integral Feynman-Kleinert formalism and the density functional Levy's one recommends the use of the Feynman-Kleinert density/densities for being implemented in density functionals with chemical relevance, the electronegativity for instance. In this regard, the density functional for Mulliken electronegativity is in next reviewed and exemplified for atomic scales.

Mulliken Density Functional Electronegativity
For an N-electronic system placed into an external potential ) (x V the general (first order) equation of the change in electronegativity, , can be written as [3,7,9,[100][101][102] in which the variation of the electronegativity χ (or the negative chemical potential in the Parr definition χ μ − = [103]) for an electronic state correlates with the number of electrons and potential variation through the chemical hardness (η ) (264), and the Fukui function (f), [3,5] being x the position vector, a.k.a. the quantum path or its average as the Feynman centroid.
Next, let us re-express the hardness and Fukui function through the relations [3,104] where S and s(x) represent the global and the local softness defined respectively as while being connected by the integral relation in turn relaying on the assumed N-normalized density, see Equation (2). Now, with Equations (385) and (386), we can integrate the equation (383) for the electronegativity to obtain under the initial natural zero electronegativity value as V(x)→0. The integrals in (390) can be carried out once the local and global softness s(x) and S, respectively, are analytically known. This can be achieved by assuming the quasi independent-particle model within density functional theory providing the softness kernel expression [104] The softness kernel (391) allows expressing the local softness s(x) by the direct integration with the effect of transforming the bi-local into the local density behavior where the well-known delta-Dirac integration rule (225b) and the density normalization condition (2) were used. Consequently, the global softness S of (389) can be immediately analytically expressed by further integrating of the local softness (392), with the result where the condition (2) accounted once more.
By introducing local and global softness expressions (392) and (393) into the Equation (390) the so called absolute electronegativity can be analytically solved yielding the formulation ( ) in terms of the additionally introduced response density functionals However, since electronegativity (394) may be identified with Parr's differential electronegativity which, at its turn can be related with the correspondent chemical Mulliken one by means of the finite difference successive transformations ( ) ( ) in terms of ionization potential (IP) and electronic affinity (EA); Yet, the integral generalization of (397) may be revealed by performing the reverse writing [7] ( ) ( ) Now, according with the Hohenberg-Kohn first theorem [1] and of the chemical action principle [7,9,18,20] the two terms in the right hand side bracket of equation (398) identically vanish since does not optimize the associations of the electronic density of one state with the external potential applied on that state, thus leaving with the identity [102] ( ) Upon the insertion of the absolute electronegativity (394) in (399) it provides the chemical Mulliken density functional [7,100,101] ( ) The exposed density functional formulation of electronegativity features reach physical contents, since the derivation appeals on fundamental quantum principles as the Hohenberg-Kohn and the chemical action theorems, while complementing somehow at the valence ( 0 → β ) level the previously density matrix one (263) -worked out in the context of (the fourth order) semiclassical expansion (271); nevertheless, this valence character of the chemical Mulliken electronegativity will be in the following tested for atomic scale through Feynman-Kleinert density implementation.
As a note, it should be mentioned that Yang has shown [105,106] how the integral formulation of the Kohn-Sham density functional theory arrives to the electronic density expression performing the Wigner semiclassical expansion combined with the short time approximation, i.e., the valence approximation, regarding the β parameter. For these reasons the present density functional application the Wigner or semiclassical limit of Feynman-Kleinert path integral approach will be used.

Atomic Electronegativity by Feynman Centroid Path Integral
The electrons of the atomic system are distinguished as the core-and the valence-ones within pseudopotential theory of atoms and molecules [100,101,107,108]; this, because it aims to provide a "valence-only" theory for these systems, while assuring the simplification of the computations. Certainly, the all-electrons picture is also possible through facing with the serious technical problem to assure the orthogonality constrains among all wave functions of all electrons of an atom. Moreover, having the valence shell treated separately is relevant for computing electronegativity, because of its definition regarding the added electron to the valence shell under the core influence. Therefore, a wise step is provided by the transformation of a many-valence electronic problem into an one-valence electronic system, so that the canonical density formulations can be at once considered. For achieving this, the link between the exact and density dependent pseudopotential is enforced by the latter's satisfying the virial theorem releasing with the radial scaling of the pseudo-orbital [108] ( ) ( ) ( ) with the scaling factor q. Therefore, the scaling factor q is searched in relation with the number of valence electrons, but such to fulfill the normalization condition  [110,111].
Within the pseudopotential methods we arrive at two possibilities for the electronic density and, consequently, for the electronegativity evaluations.
o The first one considers only the pseudo-potentials into the path integral formalism that gives the electronic density in the quantum statistical manner as it was described in the previous Section 5.4.1. This way, a strong physical meaning is assured because all the information about the electronic density and electronegativity are comprised (and dictated) only by the pseudopotential. Yet, the problem that arises in this approach is that the electronic density depends on the β parameter. This parameter will be fixed so that the electronic density to fulfill the path integral normalization condition. Additionally, the search of the β parameter must be done in the semiclassical (high temperature) limit ( 0 → β ) for which the path integral formalism corresponds to the excited (valence) states of atoms. o The second approach takes beyond to the pseudopotential data also the valence basis and the electronic densities are then computed in the accustomed quantum manner. At this point we need to consider the working orbital type for the atomic systems and we will chose the s-basis set because its spherical symmetry.
Accordingly, it follows that both electronic density approaches have their own parametric dependency. This implies that also the computed electronegativity will feature the scaling effect on the electronic density raised due to the one effective valence electronic approach. With this assumption at the background of density computation we should recover in the provided electronegativity the real (many) electronic valence state by an adequate nomination of the specific values for the β and q parameters.
At this point we need a criterion in order to properly control the re-scaling procedure. In order to unveil this criterion we are looking back on differential electronegativity formula (394) that should be seen as the kernel function for the Mulliken electronegativity functional (400). If we observe the analytical places the introduced chemical response indices a, b and the chemical action index C A appear, respectively, it can be easily seen that only the chemical action is coupled with the total number of electrons in the concerned state.
In this respect, while noting the above scaling condition (401a) as being quite restrictive for the scaling factor q, an additional constrain that takes into consideration the number of valence electrons is to be regarded. This aim is to be accomplished observing that the previous atomic electronegativity semiclassical (valence) formulation (263) fits with the definition of chemical action of Equation (5) for the Coulombic potential [115] { } with Z being the nuclear charge.
The results for the chemical action (403) computed by the two above-mentioned quantum computational schemes are collected in Table 3 with those for the electronegativity (400) in Table 4, among other significant scales. All data are comparatively represented in Figures 7 and 8.     Table 4 by path integral (PI) and basis set (BS) methods.
Analyzing these results, it is clear that for the path integral approach better correlation between the electronegativity and the chemical action trends is obtained as comparing with those arising from the s-basis set implementation. The highest discrepancy between the chemical action and the associated electronegativity within the s-basis set computation appears mostly for the first transitional row, see Figure 7.
However, we can not exclude the s-basis set electronegativity scale just through comparison between scales since other similar discrepancies appear (even in the main groups) when the Mulliken-Jaffe and the Xα methods are compared, for instance (see the Table 4). In any case, the present s-basis set results may help in judging also the various criteria of validity for an electronegativity scale, see Section 4.5.
However, it remains to show that the employed chemical action criteria (403) do not enter in conflict with the type of the orbital choice, when this is properly done. For instance, we consider the atomic systems of C, N, O with the s-and p-orbital type basis set and also the sp, sp 2 and sp 3 hybridization states. Then, by applying the re-scaling procedure according with the chemical actionelectronegativity rule (403) we get the respective electronegativities and chemical actions for both the path integral and basis set implementations using the pseudopotential data [109], as shown in Table 5 and drawn in Figure 9.   Table 5.
They, nevertheless reveal how close the values of the chemical actions and corresponding orbital electronegativities are, in general. Such feature is susceptible for extension in treating the chemical bonds as was recently employed [20,116].

Levels of Non-equilibrium Dynamics
The electronic states associated to the cyclic reactions, oscillatory chemical phenomena, systems with instabilities, etc., are modeled by the dynamics of discrete states subjected to the general master equation [117][118][119][120][121] or, for the continuous transformation of states -by the quantum evolution equation where W or W n represent the discrete or continuous probability density, respectively, while w describes the transition probability; yet, when it is given by the (diffusion) distribution the filtration property of Dirac functions (225b) and the obvious relationships In essence, Fokker-Planck equation reflects the evolution of the quantum fluctuations (or of the probability density) under the external applied potential ) (x V by individuating the drift factor while being accompanied by the diffusion factor ) ( also known as the stochastic noise [123][124][125].
Such description of open electronic systems corresponds with the probability density level of nonequilibrium evolution. It originates in modeling of the so called Markovian processes by the reactive chain with correlated successively events. The pairing of the neighboring events is accommodated within the density matrix, the quantum propagator, or the conditioned probability density ( ) It is worth noting that for the Markovian processes all relevant information is contained in the first two probability functions 2 1 &W W being these involved in all successive correlated events; as such, with the properties (410) they may produce any desired extended propagator. For example, one could be immediately write that which, while equated with the equivalent one in close correspondence with the group property of propagators (141). Now, in terms of conditioned probability density, i.e., by using the property (410c), the Fokker-Planck equation (407a) rewrites as Very interesting, the Fokker-Planck equation (413) may be resumed under the hydrodynamic form through introducing the current of probability density On the other hand, for the initial condition of probability density identified with the Dirac distribution the property (410c) yields the predicted probability density as the working propagator form. Alternatively, if the initial condition is somehow relaxed to the Gaussian form to which any continuous distribution is reduced according with the central limit theorem, then, the same general property (410c) prescribes the Markovian probability density of the type as a sort of smearing for the quantum propagator, being this behavior the specific characteristic for the non-equilibrium systems where fluctuations should be averaged in order to produce predicted probabilities. In either of cases, Equations (420b) or (421b), the working Markovian probability current (415) unfolds like staying as the basis for the analytical representation of the non-equilibrium electro-reactive dynamics in open environment, as atoms-in-molecules are, for instance.

Non-equilibrium Lagrangean
Once identified with the quantum propagator, the above probability density may be readily expressed by a special form of the Feynman path integral; namely it may be viewed as the quantum stochastic Fokker-Planck (FP) propagator [126,127] marking the appearance of the third form of path integral representation, after those specific to the quantum mechanics and quantum statistics, see Equations (121a) and (121b), respectively. Nevertheless, the stochastic Lagrangean in (423) assumes the general ansatz [48] ( ) ( ) ( ) Now, employing the (425b) expansion the potential dependencies in (426) may be as well considered until the second order cut-off series so that the probability density property (410c), forgetting the subscript -since no reason for confusion, becomes: where we used the notations Next, based on the observation that it follows the coordinate displacement proportionality which, along the temporal abstracted counterpart from Equation (425a) allows the first order expansion in ε for the elementary probability density With these rules, the expression (433) simplifies to or even more as Finally, worth remarking that because the quantum Planck constant h is apparently absent in the Fokker-Planck approach -it seems more appropriately called as mesoscopic picture, thus laying between the quantum microscopic and Newtonian macroscopic pictures of natural systems. It nevertheless works for assessing new insight of the electronic behavior at the atomic and molecular levels by means of the averaged harmonic and anharmonic fluctuations -as is to be revealed below.

Harmonic Markovian Density Matrix
Harmonic solution of non-equilibrium propagator is to be found by using the method of path integrals, starting from the associated Lagrangean explicated in the previous paragraph. It starts from the expression of conditioned density of probability in the form of path integral (423) with (439) [48] ( ) that may be formally simplified by considering the spectral shift of eigen-FP values introducing the potential transformation The form (442) is next employed for the Onsanger-Matchlup harmonic potential (equivalent with the Brownian motion) [ obtained from equation (413) with the drift Equation (407b) re-shaped with the potential (443), while considering for the diffusion coefficient ) ( the constant D.

Anharmonic Markovian Density Matrix
The next challenge stays the finding of non-equilibrium propagator for the anharmonic generalized potential Viewed as a perturbation, the transition probability or the stochastic propagator (442) will be calculated by neglecting all the contributions 2 , ≥ n g n , successively as [48] ( ) The involved classical (235) and Green function (248) propagators are rewritten through the correspondences (445) with the factor m / h of (233b) toward the stochastic picture the Heaviside function (247). Now, while employing the Wick correlation rules (220) and (221) for the identical end-point times, one gets the Feynman expressions as combinations of the components of Equations (453).
with the coefficients (456a). The last cross-check regards the stationary solution behavior. This is done by verifying the identity between the stationary anharmonic solution abstracted from Equations (455) and (456), which is With the help of frontier (Dirichlet) constrains have been employed with the forms arisen from combining the relationships (407b) and (441) upon it followed by its replacement in the solution (463a) the identical form with that of Equation (460) is provided, thus certifying the correctness of the present general solution (455)-(456) also in the stationary (i.e., the asymptotic limit) of anharmonic potential (449).
Remarkably, the asymptotic limit of the anharmonic solution of the Fokker-Planck case displays a form without the memory of the initial event ( ) a a t x , , thus susceptible to be assumed as a "special" canonical density for specific electronic analysis. In which way? It will be unfolded in the next section.

From Thom's Catastrophe Concepts to Chemical Bond Topology
The Thom's catastrophe theory [129] basically describes how, for a given system, a continuous action on the control space (C k ), parameterized by c k 's, provides a suddenly change on its behavior space (I m ), described by variables x m 's, through the stable singularities of the smooth map being η(c k , x m ) called the generic potential of the system. Therefore, catastrophes are given by the set of critical points (c k , x m ) for which the field gradient of the generic potential vanishes or, more rigorously: a catastrophe is a singularity of the map M k×m → C k .
Next, depending on the number of parameters of space C k (named also as the co-dimension, k) and of the number of variables of space I m (named also as the co-rank, m), René Thom had classified the generic potentials (or maps) given by Equation (465) as seven unfold elementary (in the sense of universally) catastrophes, i.e., providing the many-variable (with the co-rank up to two) -manyparametrical (with the co-dimension up to four) polynomials, listed in the Table 6. Going to the higher derivatives of the generic potential (the fields), it will be said that the control parameter c k * for which the Laplacian of the generic potential vanishes gives the bifurcation point. Consequently, the set of control parameters c # for which the Laplacian of a critical point is non-zero defines the domain of stability of the critical point. There is clear now that the small perturbations of η (c*, x) bring the system from a domain of stability to another; otherwise, the system is located within a domain of structural stability. Remarkably, the above described cases correspond to the equilibrium limit of a dynamical (nonequilibrium) evolution for an open system 0 ,... ) ; ( ); ; where the behavior space is further parameterized by the temporal paths x m (c k , t). The connection with equilibrium is recovered through the stationary time regime imposed on the critical points. This way, the set of points giving a critical point in the stationary +∞ → t regime (the so called ω-limit) corresponds to an attractor, and forms its basin, whereas the stationary regime −∞ → t (the so called α-limit) describes a repellor.
These catastrophe concepts have the merit to describe the evolution of local properties of (in principle) any natural system. In this framework, the chemical bonds can be seen as the equilibrium part of the evolutionary binding processes. Therefore, to describe the bonds and binding, a suitable generic potential η(x) has to be consider. In topological studies of electron localization across a chemical reaction modeled by the catastrophe approach the variable x can stay also as the reaction coordinate.
With the aim of properly choosing the function η(x), with x the space-spin coordinate, the electronic wave function Ψ, as provided by Schrödinger or related Hartree-Fock formalisms, can be an option, but suffers from the lack in the real space significance. Next, a better choice regards the electronic density ρ(x) as the real descriptor of the topological electronic distribution in space and for identifying the bonds as well. In this respect, the gradient equation of the electronic density, ∇ρ = 0, provides the critical points whereas the Laplacian equation, Δρ = 0, indicates the bifurcations and stability zones, respectively. This picture was intensively used by the Bader's theory of atoms in molecules [84], with partially success. An alternative electronic topological approach was performed by Mezey by considering the changes in shape of the bonding isosurfaces [130], but in a form that does not allow the description of bonding in terms of Laplacian. Worth noting here that the Laplacian plays a crucial role in topological bond description, being associated with the quantum mechanical transcription of the kinetic electronic energy, h , being at its turn related with the minus of total energy of the electronic system, E = -T, through the virial theorem at equilibrium [108].
In fact, this feature of Laplacian was extensively employed by Bader's atoms in molecules theory at the purely electronic density level in order to quantum rationalize the previous Gillespie's geometrical VSEPR (Valence Shell Electron Pair Repulsion) description of the molecular bonds [131,132]. A more elaborated choice in generic potential was proposed by Becke and Edgecombe through introducing the electron localized function (ELF), representing a density combination rather than the electronic density solely. Their approach (abbreviated as "BE") prescribes η(x) with the form [49] and: as a combination between the kinetic electronic density terms from the Hartree-Fock (or Kohn-Sham) orbitals φ i , the Weizsäcker gradient correction, and the homogeneous (abbreviated by "h") Thomas-Fermi descriptions, respectively [133]. In Equation (469) D(x) accounts for the excess of local kinetic energy density due to Pauli repulsion, whereas D h (x) plays the role of the "renormalization" factor. However, this ELF function behaves like a density by mapping its values onto the realm [0,1], where 1 corresponds to the perfect electronic localization, being therefore suitable for the gradient and Laplacian performances. The recent topological studies have revealed the effectiveness of the above η BE (x) function in describing both the electronic localization in bonding as well for modeling the chemical reaction pathways [134].
However, despite of η BE (x) efficiency in bonding characterization, a series of aspects regarding its appearance in the context of the universal unfolding of the catastrophes (see Table 6) as well as within the time-dependent and the stationary ω-limit have remained unexplored. The present Markovian description of the anharmonic potentials with the help of Fokker-Planck equation and of its path integral solution aims filling this gap -as exposed next.

Fokker-Planck Approach of Electron Localization
In order to better understand the actual ELF approach, it is worth reminding that the origin of the above η BE (x) relies in evaluation of the conditioned pair probability with which one electron is located at point x b with the spin σ once the reference electron is located at point x a with the same (parallel) spin σ [135] ( ) being the coefficient σσ A identified with the function D(x) in Equation (470) within the so called "hole" function approach, whereas in the present treatment has to be re-determined. Nevertheless, η BE ELF emphasizes on the key role of the conditioned probability density -to be here considerate alternative Markovian description of natural processes. As previously shown, the Markovian treatment of the conditioned probability density given by the general Fokker-Planck path integral (442) with correspondences (445) in atomic units see its definition from (437) and the transformation (441). Within the anharmonic potential case it features the non-linear shape assuring the connection with the electronic localization by the homogeneous and inhomogeneous (or gradient) specializations with the help of electronic functions (471) and (470), respectively. Yet, the correspondences (475) are motivated as follows [50]. If one considers the working effective potential (449) as the bilocal dependency it models the field produced by the reference electron located at x a over its spherical neighborhood which contains the coupled electron at the distance (or radius) x b , thus being characterized by the density (radial) equation of Poisson type while clearly revealing the role of the homogeneous and gradient related terms as being the friction γ and perturbation factor g in (476), respectively. Note that the form of the potential (476) assumes one such that the new parameter g* to be strictly positive, while being a transformation of the negative of g, along the general homogeneous h parameter dependence. Finally, with Equations (487) and (488) back in the founded bifurcation parameters A σσ , and further into the unfolded catastrophe functions (483) and (484), the corresponding general Markovian ELF cases are displayed as [50]  i.e. displaying an inverse function of a gradient to homogeneous ratio electronic contributions (470) and (471), respectively. Such general definition immediately allows its physical interpretation. As such, when density gradient dominates, i.e., Resuming, the right definition of ELF has to provide the limits [50] corresponding to the minimum/maximum error in spatial localization, respectively, and for which the interpretation has to indicate where the electrons are trapped rather than where they have peaks of spatial density, as is commonly interpreted [5,15]. In short, the meaning of ELF is associated with the error in spatial localization of electrons, being zero when the electrons are precisely located, through observing their density gradient distribution.

Working Markovian ELFs
For practical purposes the general Markovian ELFs given by Equations (490) Nevertheless, the ELF cases given by Equations (495) and (496) correct, in a Markovian framework, the previously purely hole-pair probability approach. The question is to decide which of the above Markovian ELF' cases are more "corrective" respecting the Becke-Edgecombe one. For better visualizing the answer the Figure 11 shows the BE-M1 and BE-M2 differences on the relevant homogenous (h parameter) versus inhomogeneous (g parameter) contributions to electronic localization.
The analysis of Figure 11 clearly reveals that for a moderate inhomogeneous contribution to the electronic gas the first Markovian ELF of Equation (495) corrects the Becke-Edgecombe localization function up to 15%, whereas, in the same conditions, the second Markovian ELF of Equation (946) improves only up to 8% the Becke-Edgecombe ELF treatment. Therefore we can conclude that the most corrective Markovian ELF to the Becke-Edgecombe approach stays the first dependence of Equation (495); this one can be further tested for prediction of the electronic localization in atoms and of bindings in molecules. They produce the corrections up to 30% and 20% for the Becke-Edgecombe ELF approach, respectively, in the moderate inhomogeneous electronic behavior -as the Figure 12 reveals. Again, the Markovian ELF corresponding to the first case, Equation (498), is the most corrective respecting the Becke-Edgecombe treatment. The last considered ELF particularization eventually involves the hyperbolic trigonometric function with the form: x D Now, the differences to localization introduced by these Markovian ELFs as referring to the Becke-Edgecombe formulation are analyzed through the representations given in Figure 13; the analysis sharply indicates that the Markovian ELF++ approaches depart between 10-20% from the Becke-Edgecombe ELF, providing an intermediary situation between Markovian ELF (8-15%) and Markovian ELF+ (20-30%) predicted by Equations (495)-(496) and (498)-(499), with the representations of the Figures 11 and 12, respectively.
Overall, judging by both the analytical complexity and meaningful physical background, grounded on the Fokker-Planck approach of the non-equilibrium towards equilibrium systems, we propose that the Markovian ELF1+ of Equation (498) to be adopted as the electronic localization function for the practical topological characterization of the atomic shells and the molecular bonds [136][137][138][139][140].
The combined path integral with the non-linear and electronic density aspects fully qualify our analytical results as a reliable framework within which the electronic localization targeting the bonding evolution theory to be further developed in the years to come.

Conclusions
Since the recent most celebrated quantum theory of Chemistry -the Density Functional Theory (DFT) is mainly based on density functionals, which relay on their turn on the many-body densities, the search for electronic density both as computational (here understood as analytical) and conceptual assignments remains a crucial endeavor in quantum chemistry, comparable with the landmark theoretical predictions of spectra in the early age of quantum physics.
Yet, for achieving such challenging task the complex mathematical-informatics and mathematicalphysics seems to be at the foremost background for computational and conceptual Density Functional Chemistry, respectively. The present review was dedicated to the latter goal that is to present the analytical framework in which the many-electronic systems may be described by the associate densities at various levels of conceptualization, approximation, and applications.
As such, through presenting the basics of density matrix theory, the precursor of DFT, the path integral concept appears as the natural solution for expressing the time-space electronic density. Indeed, the Feynman path integral formulation has been revealed as the natural generalization of the Schrödinger equation, being in close relation with the propagators and Green function of a given quantum system, either at equilibrium or coupled with a temperature bath or particle environment.
Nevertheless, the density matrix -path integral description allows the general formulation for the many-electronic density through the so called canonical density algorithm; it prescribes that the system is firstly solved for the single electron evolution under the concerned potential for which the timespace density matrix is analytically formulated, in an evolution manner, as the propagator Remarkably, the partition function involvement in this density algorithm was widely and most extensive used by the Feynman-Kleinert approach which was proved to furnish meaningful approximations either for the ground state (as was the case for atomic Hydrogen and the Bohr's orbital proofed stability) as well for the higher temperature or excited or the valence states (that resembles the semiclassical approximation).
Regarding the realization forms of path integral approaches of a quantum problem/system there were individuated three major pictures: the quantum mechanical (QM), the quantum statistical (QS) and the Fokker-Planck (FP) ones; the passage among them as well as their inter-conversion and equivalence being realized through the time transformations presented in the Table 7. It, nevertheless, leads both with philosophical and practical consequences: epistemologically, there seems that the time itself may suffer transformations being of quantum order ( h ) and correlated with inverse of the thermic energy ( β ), while passing from quantum to statistical description of Nature; as well, the mass in equilibrium states plays the role in open systems of inverse of diffusion (naturally, since presenting inertia) but also with a quantum manifestly nature, while the ordinary (QM or QS) harmonic oscillations' frequency becomes friction in Fokker-Planck description of non-equilibrium systems.
Moreover, practically, such inter-conversion table allows for immediately transferring of one result obtained within a quantum picture into another without the need of entirely problem reformulation. This procedure was largely considered and applied throughout the present review; it nevertheless leads with another epistemological conclusion, namely that the quantum mechanical oscillatory description is equivalently converted into the hyperbolic function for statistical and Markovian (or FP) frameworks, which further means the quantum modeling by the Gaussian wave-packet and Green functions. At this point, worth being mentioned the explicit proof for the de Broglie equivalence with the Gaussian wave function by the smearing out procedure of the fluctuation of the closed paths with the effective partition function approach; even more, such equivalence is justified by the very roots of quantum theory since the Born normalization of the de Broglie wave-packet is finely satisfied by the Gaussian form of its Fourier coefficients (amplitude) [67].
From the chemical point of view, the valence states are those situated in the "chemical zone"-and they are the main concern for the chemical reactivity by employing the frontier or the outer electrons; consequently, the semiclassical approximation that models the excited states was expressly presented either as an extension of the quantum Feynman path integral or as a specialization of the Feynman-Kleinert formalism for higher temperature treatment of quantum systems. However, due to the correspondences of Table 7 one may systematically characterize the semiclassical (or quantum chemical) approaches as one of the limiting situations: : the quantum statistical short-time limit; o ∞ → T : the high-temperature limit; o 0 → ω : the flat potential, or the quasi-homogeneous (Thomas-Fermi) limit; yet, this may be easier visualized by noting the discrete-to-quasi continuum transformation of eigenlevels intervals in the exited zones of quantum systems (atoms, molecules), i.e., where the approximation 1 ) /( << = T k B h h ω β ω holds; moreover, this limit nicely overlaps with the "free harmonic approximation" used in this work, when the interplay between the free and harmonic motion helped in elucidating and solving (by integrating out) the quantum fluctuations along the classical paths; : the bosonic limit due to the scaled equivalence 3 / 1 N T or 3 / 1 Z T when the system is thermally expanded, being it related with the Thomas-Fermi theory, here not exposed. This is the way the chemical reactivity basic indices of electronegativity and hardness were computed both by a direct semiclassical density matrix implementation as well as by the author's derived density functional of Mulliken electronegativity within Feynman-Kleinert specialization to higher temperatures; the results had confirmed that the direct density matrix formulation for chemical reactivity indices provides almost not-observable values for the valence states with high quantum number, although nicely ordering the periodicity trends across periods of elements, while the more elaborated density functional pseudo-potential path integral & basis set implementations into the actual density functional Mulliken electronegativity produces much closer values to those computed upon the observed ionization potentials and electronic affinities. Finally, the path integral formulation of the Markovian open systems leads with general formulations for the electronic localization functions, enlarging this way the implementation possibilities beyond that previously derived by Becke and Edgecombe for isolated (at equilibrium) quantum systems.
It is therefore hoped that future cross-fertilization between path integrals' physical approaches and the chemical concepts of reactivity will conduct towards the unification of the chemical bonding types and phenomena within a theory of the chemical field. and promoters of Path Integral method worldwide, to whom I am indebted for many shining ideas on path integral formulations during my last decade repeated visits to the Free University Berlin and the University of Duisburg-Essen, besides various meetings with occasion of international path integral and physical-mathematical conferences, and by their most recent honoring, inspiring and stimulating visit to West University of Timisoara. Finally, this work is a humble homage to the Feynman spirit and to what it helps to modeling in Natural Sciences.

Appendix: Poisson Formula for Series
Aiming the evaluation of series