Modern State of the Conventional DFT Method Studies and the Limits Following from the Quantum State of the System and Its Total Spin

Abstract

At present, the density functional theory (DFT) approach became the most widely used method for study molecules and solids. In the atmosphere of such great popularity, it is particularly important to know the limits of the applicability of DFT methods. In this chapter, I will discuss the modern state of DFT studies basing on the last publications and will consider in detail two cases when the conventional DFT approaches, in which used only electron density and its modifications by gradients, cannot be applied. First, the case related to the total spin S of the state. As I rigorously proved for an arbitrary N-electron state by group theoretical methods, the electron density does not depend on the total spin S of the state. From this follows that the Kohn-Sham equations have the same form for states with different S. The critical survey of elaborated DFT procedures, in which the spin is taken into account, shows that they modified only exchange functionals, and the correlation functionals do not correspond to the spin of the state. The point is that the conception of spin in principle cannot be defined in the framework of the electron density formalism, and this is the main reason of the problems arising in the study by DFT approaches the magnetic properties of the transition metals. The possible way of resolving spin problems can be found in the two-particle reduced density matrix formulation of DFT. In the end, it will be considered the case of the degenerated states, in which, as follows from the adiabatic approximation, the electron density may not be defined, since electronic and nuclear motions cannot be separated, since, the vibronic interaction mixed them.

1. Introduction

The DFT method developed by Walter Kohn with his collaborators more than 55 years ago is widely used for study molecular systems and solids. First, Kohn published with Hohenberg [1] their famous theorem on which the DFT theory is based. Then Kohn with Sham [2] obtained their well-known Kohn-Sham equation. The number of citations on these first Kohn and collaborators publications is increased every year [3]. In 2010 papers [1, 2] were cited 11,000 times and in 2012 14,000 times. So, we can expect that at present the number of citations will be larger.

Last years, the number of Congresses on Nanosciences and Nanotechnologies is considerably increased. In one of the reports made in Las Vegas on October 2016 by Miyazaki, it was claimed that using the modern computational facilities they were able to apply the DFT molecular dynamic simulation to systems with million atoms, see Ref. [4]. Certainly, it can be done only using the traditional Kohn-Sham (KS) formalism based on the electron density $\rho \left(r\right)$ and its modifications by gradients.

The electron density is the diagonal element of the spinless one-particle reduced density matrix,

$$\rho \left({\mathit{r}}_1\right)=N\sum _{{\sigma }_1,\dots ,{\sigma }_N}\int {\left|\mathrm{\Psi }\left({\mathit{r}}_1{\sigma }_1,\dots ,{\mathit{r}}_N{\sigma }_N\right)\right|}^{2}d{V}^{\left(1\right)}$$E1

In Eq. (1), the spin projections ${\sigma }_i$ span over the whole spin space and the spatial coordinates are integrated over the N − 1 electrons excepting the first. If Ψ is defined in 4N-dimensional Gilbert space, $\rho \left(r\right)$ is defined in a three-dimensional space. Evidently, the calculations, in which only $\rho \left(r\right)$-formalism is used, will be considerably faster than Ψ-formalism is used, and the $\rho \left(r\right)$-formalism can be applied to larger systems.

Let us consider another point: what we lose, if we use ρ-formalism with the modulus of ${\left|\mathrm{\Psi }\right|}^2$ squared, as in Eq. (1), instead of the wave function formalism. It is evident that in the transition from the wave function formalism to the probability density ${\left|\mathrm{\Psi }\right|}^2$, we lose the phase of the wave function. Due to the insensitivity of the probability density to the symmetry of the state (we will discuss it in the following sections), we also lose the symmetry characteristic of the wave function and cannot determine the Pauli permitted states, on which molecular spectroscopy is based. The diagonal element of the full and all reduced density matrices, as it was proved in my studies and discussed in book [5], does not depend on the symmetry of the state and its dimensionality.

It can be expected that after integration we lose some information. In the case of electron system, the one-particle reduced density matrix must be used. This leads to loss of information connected with the two-particle correlations, which are described by the two-particle reduced density matrix.

What is not evident and deserves a special discussion: as it was proved by the author for an arbitrary many-electron system, the total spin S of the system in principle cannot be introduced in the DFT studies. This can be done at the two-particle reduced density matrix level. In the Section 2, we will discuss this problem in detail.

On the other hand, even at the framework of the two-particle reduced density matrix formalism, one cannot study the non-additive many-body effects, which determined by many-body forces, In this connection, I would like to mention that when I arrived from Moscow to Mexico by invitation of the Director of IF-UNAM Octavio Novaro and was working in his laboratory, we obtained a closed formula for the energy of N-body interactions [6].

Later on, in the author book [7], Chapter 4, different general cases for the many-body forces have been considered, see also paper [8]. In several publications, e.g., in Refs. [9, 10], it was conclusively demonstrated that the clusters built from closed-shell atoms (atoms without valence electrons) are stabilized by the three-body forces. The alkaline-earth clusters: Be_n, Mg_n, Ca_n, and etc. are the typical example of such clusters. They are stabilized by the three-body forces, It is also important to mention that for the stability of rare-gas clusters the three-body dispersion forces, which are known as the Axilrod-Teller-Muto forces, play a decisive role, see Section 4.3.3 in Ref. [7] and recent review by Johnson and co-authors [11].

It is instructive to discuss shortly the physical sense of the non-additivity and many-body forces concepts. As is well known, the interaction of charge particles is described by the Coulomb law

$$V=\sum _{a<b}\frac{q_a{q}_b}{r_{\mathit{ab}}},$$E2

where charges $q_i$ are considered as points. Eq. (2) contains only two-particle interactions, so, it is additive. However, in quantum mechanics, the charge particles are not points and they are not rigid. Atoms and molecules obey quantum-mechanical laws. If we consider the Coulomb interactions between charged atoms (or charged molecules), the additivity is lost. The interaction will depend on surrounding. In the case of three atoms, the third atom can polarize the electronic structure of two others, and this leads to three-body forces, since the interaction depends on three interatomic distances.

It should be noted that the possibility of application of DFT approaches to large systems, which were not available to be studied before, induced a euphoria in the DFT community. This euphoria led to wide using DFT methods without an analysis of the limitations following from quantum mechanics.

In many publications it was revealed that the applications of DFT method in some cases lead to incorrect results. First, it was recognized in the DFT studies of intermolecular interactions. The potential curves obtained by the early created DFT functionals for many stable in experiments dimers were repulsive, since in these DFT functionals the dispersion energy was not taken into account.

Then it became clear that DFT methods meet serious difficulties in studies of transition metals with nd electrons. These problems were analyzed by many authors, e.g., by Cramer and Truhlar [12]. I would like to stress that most of difficulties discussed in their review [12] are connected with the problem of spin in DFT approach and in principle cannot be resolved in the framework of the electron density, ρ(r), which belong to the one-particle reduced density matrix, see Eq. (1). The spin problem in DFT approaches will be analyzed in detail in Sections 2 and 3.

Last years, many comparative studies of the relative precision of exchange-correlation (XC) functionals are published. Below I will discuss some important, from my viewpoint, DFT papers published on this topic in the last years. Certainly, the list of selected papers is only a little part of thousands DFT papers that are published each year.

Gillan et al. [13] analyzed different kind of XC functionals for liquid and ice water and water clusters. The conclusion was that many functionals are not satisfactory because do not describe correctly the dispersion. Let us stress that this situation takes place for such widely studied substance as water. The authors mentioned that after they included in XC functionals the non-local dispersion, the results still cannot be admitted as completely satisfactory.

In the Taylor et al. [14] paper, the precision of DFT calculations of intermolecular interactions with respect to highly accurate benchmarks for 10 dimers was analyzed. Their review is comprised 17 authors, among them are the well-known creators of XC functionals Angyán, Hirao, Scuseria, Truhlar and others.

The perspectives of DFT theory were discussed in the paper by Truhlar and collaborators [15]. The authors also analyzed recent Minnesota functionals. On the other hand, Mardirossian and Head-Gordon [16] benchmarked the Minnesota functionals using a very comprehensive database and came to conclusion that none of them are state-of-the-art for non-covalent interactions and isomerization energy.

In the article published by Medvedev et al. in collaboration with Perdue [17], 128 XC functionals created in period 1974–2015 were analyzed. The authors [17] made the comparison of the normalized error for different functionals for atoms and its ions, see Figure 1, using the CCSD approximation as a reference level. It is important to stress that according to their results the normalized error of $\rho \left(r\right)$ up to 2000 decreased and then increased due to the introduction of semiempirical functionals.

Figure 1.

The precision of calculating ρ using different functionals.

Kepp [18] made a comment on the conclusion by the authors [17] that after the early 2000 the functionals strayed from the path toward exactness. Kepp indicated that the “straying” is not chemically relevant to the studied systems. In their response [19], the authors performed a special calculation, from which follows that their results are valid also for heavier system, including molecules.

In the next paper published by Perdew with Medvedev et al. [20], the authors discuss the possibilities and pitfalls of statistical error analysis, taking as an example the ranking of approximate functionals by the accuracy of their self-consistent electron densities.

As mentioned above in our discussion of last publications, the DFT community is concerned on the still existing problems in some applications of DFT approaches and most of these problems are connected with the quality of XC functionals. I like and completely agree with the witty comment made by Savin [21]:

“The existing great number of different expressions for the XC functionals can be considered as evidence that we still have not satisfactory correct XC functionals”.

Among new publications, it is important to mention the very recent paper published by Perdue with collaborates [22], in which the problem of the symmetry breaking in DFT was discussed. The interesting paper was recently published by Bartlett [23]. Even the title of the paper “Adventures in DFT by a wavefunction theorist” looks quite intriguing. I also recommend readers the quite comprehensive and interesting review by Johnson and co-authors [11] published in 2021.

In my following discussions, I will analyze in detail the two problems:

The first problem is connected with the definition of spin in the KS-DFT framework. It will be shown that the concept of the total spin S of the state, in principle, cannot be defined in the frame of electron density formalism.

The second problem is related to the degenerate states in DFT. As follows from the general conceptions of quantum mechanics, at Born-Oppenheimer approximation in the case of degenerate states, the electronic and nuclear motions cannot be separated, they are mixt by so-called vibronic interactions. This problem will be discussed in Section 3.

In this chapter I will not discuss Ψ-versions of DFT that solves many problems in DFT but lost the simplicity of KS-DFT formulation. The Ψ-formalism was successfully combined with DFT approach by Gőrling, Trickey, and some other investigators. Nevertheless, the existing problems in the applications of DFT approach are still remaining actual. I will discuss them in the next two sections. The discussion will be based on two theorems proved by the author in 2007 [24].

2. Electron density ρ and the total spin S

It was recognized sufficiently long ago that the concept of the total spin of the many-electron quantum state is quite difficult to base in the frame of DFT approaches. Studying the two-electron system, McWeeny [25] came to conclusion the electron density does not allow to identify the spin state. McWeeny formulated it as the following statement:

“Electron spin is in a certain sense extraneous to the DFT”.

In their analysis of DFT foundation, Weiner and Trickey [26] came to conclusion that

“…the way that the ρ-based XC potential takes account of spin is very obscure except in the simplest configurations”.

The statements of McWeeny and Weiner-Trickey are quite cautious. I would like to formulate it more definite:

“The conception of spin in principle cannot be defined in DFT at the level of the first reduced density matrix”.

To the best of my knowledge, the concept of spin was discussed in DFT community only for two-electron systems. In my article [24], I proved the theorem that the electron density does not depend upon the total spin of the state for N-electron system. This proof was done applying the formalism of the permutation group.

I would like to mention that in the $\mathrm{\Psi }$-formalism used by Slater [27], he presented the wave function as a linear combination of determinants, corresponding to a given value of the total spin S. It does not allow to obtain any conclusions about spin. On the other hand, as I noted above, in my articles, firstly in Ref. [24], I obtained the proof for N-particle system for the independency of the electron density from the total spin S of the state. It was achieved applying the permutation group apparatus. For understanding the following text, I presented in this section the short description of the formalism of permutation groups. For more detail description, see Appendix in the end of this chapter or chapter 2 in my book [28].

The permutation group is characterized by Young diagrams $\left[\lambda \right]$:

$$\left[\lambda \right]=\left[{\lambda }_1{\lambda }_2\dots {\lambda }_k\right],$$

$${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_k,\sum _{i=1}^k{\lambda }_i=N.$$E3

The rows in the Young diagram are represented in decreasing order. The presence of several rows with equal length ${\lambda }_i$ is indicated by a power of ${\lambda }_i$. For example, [λ] = [2²1²] is depicted graphically as

I would like to mention that in Appendix, I describe the unusual biography of Alfred Young who was not a professional mathematician but was a country clergyman.

Though the concept of spin has enabled to explain the nature of chemical bond, the electron spins are not involved directly in the formation of the latter. The interactions responsible for chemical bonding have a purely electrostatic nature. In non-relativistic approximation, the Hamiltonian does not depend on the spin, that means that the spin is saved, and we can operate with the value of the total spin S in the considered state.

According to PEP, the total electron wave function can be constructed as a sum of product of the spatial and spin wave functions symmetrized in respect to the irreducible representations ${Г}^{\left[\lambda \right]}$ [28].

$${\mathrm{\Psi }}^{\left[1^N\right]}=\frac{1}{\sqrt{f_{\lambda }}}\sum _r{\mathrm{\Phi }}_r^{\left[\lambda \right]}{\mathrm{\Omega }}_{\stackrel{\sim }{r}}^{\left[\stackrel{\sim }{\lambda }\right]}.$$E4

In Eq. (4), $\left[\lambda \right]$ is the Young diagram and ${Г}^{\left[\stackrel{\sim }{\lambda }\right]}$ denotes the representation conjugate to ${Г}^{\left[\lambda \right]}$. Its matrix elements are

$${\mathrm{\Gamma }}_{\stackrel{\sim }{r}\stackrel{\sim }{t}}^{\left[\stackrel{\sim }{\lambda }\right]}\left(P\right)={\left(-1\right)}^p{\mathrm{\Gamma }}_{\mathit{rt}}^{\left[\lambda \right]}\left(P\right).$$E5

where p is the parity of permutation P. The spin Young diagram $\left[\stackrel{\sim }{\lambda }\right]$ is dual to $\left[\lambda \right]$, i.e., it is obtained from the latter by replacing rows by columns. For example,

Let us return to Eq. (4), where the sum is taken over all basis functions of the representation. The normalization of the total wave function is provided by the factor $1/\sqrt{f_{\lambda }}$. It should be mentioned that the electron spin has only two projections $s_z=\pm \mathit{½}$, therefore the spin Young diagram $\left[\stackrel{\sim }{\lambda }\right]$ must have no more than two boxes per columns. In one box, the projection of spin $s_z=\mathrm{½}$ in the other box $s_z=-\mathit{½}$. It means that the total spin of this column equals 0. It is evident that the contribution to the total spin of the system of electrons will come only from uncoupled electron spins, that is, from the length of rows. The value of the total spin corresponding to spin Young diagram $\left[\stackrel{\sim }{\lambda }\right]$ is equal to

$$S=\frac{1}{2}\left({\stackrel{\sim }{\lambda }}^{\left(1\right)}-{\stackrel{\sim }{\lambda }}^{\left(2\right)}\right).$$E6

Eq. (6) enables one to find easily the values of the spin S for each spin Young diagram. For example, the spin corresponding to the spin Young diagram $\left[\stackrel{\sim }{\lambda }\right]=\left[31\right]$, is equal to $S=1$.

Let us mention that in the case of particles with $s>\mathit{½},$ for a given Young diagram can correspond several values of $S.$If spin of the particle $s=1$, to $\left[\lambda \right]=\left[31\right]$ can be attributed three values of $S=1,2,$ and $3$, see Table 2 in Section C4 in Appendix C of my book on PEP [5].

As follows from experiment, the wave function of elementary particles can be only completely symmetric or antisymmetric. It allowed to generalize the PEP, primary formulated by Pauli only for electrons, for all elementary particles:

The only possible states of a system of identical particles possessing spin s are those for which the total wave function transforms upon interchange of any two particles as

$$P_{\mathit{ij}}\mathrm{\Psi }\left(1,\dots ,i,\dots j,\dots ,N\right)={\left(-1\right)}^{2s}\mathrm{\Psi }\left(1,\dots ,i,\dots j,\dots ,N\right).$$E7

That is, it is symmetric for integer values of s (the Bose-Einstein statistics) and antisymmetric for half-integer values of s (the Fermi-Dirac statistics).

According to Ehrenfest and Oppenheimer [29], this formulation is valid not only for elementary particles, but it is valid for different composite particles as well. As examples of composite particles, the authors considered atoms, molecules, and nuclei composed by electrons and protons (at that time the neutron had not been discovered). According to the presented above the general formulation of PEP for elementary particles, the wave-functions that described them can have only two types of symmetry: completely symmetric or antisymmetric, depending on their intrinsic value of spins.

The composite particles considered by Ehrenfest and Oppenheimer [29], were composed by fermions, that is, from particles with spin ½. So, the even number of particles leads to the Bose-Einstein statistics and odd number to the Fermi-Dirac statistics,

To the best of my knowledge, the scientists that had developed methods allowing to use the conception of spin in DFT calculations considered only two-particle systems. The general case of the N-electron system was considered firstly by the author in Ref. [24] where it was proved the theorem named by some authors as the Kaplan Theorem 2. This theorem was formulated in the following manner:

“The electron density of an arbitrary N-electron system, characterized by the N-electron wave function corresponding to the total spin S and constructed on some orthonormal orbital set, does not depend upon the total spin S and always preserves the same form as it is for a single-determinantal wave function.”

According to this theorem, for any permutation symmetry of the spatial wave function described by the Young diagram [λ] that correspond to a definite value of spin S, the electron density is equal

$${\rho }_t^{\left[\lambda \right]}\left(\mathit{r}\right)=\sum _{n=1}^N{\left|{\phi }_n\right|}^2.$$E8

It is a well-known expression of the electron density for the state described by the one-determinantal function with single-occupied orbitals. It can be shown that in the case of orbital configuration with arbitrary occupation numbers, the final expression (6) has not changed and will also correspond to the electron density for the one-determinantal function.

As follows from the discussion in the beginning of this section, at the first reduced density matrix approximation the concept of spin in principle cannot be introduced in the frame of traditional KS approach and at the gradient correction level as well. In more detail see the discussion based on the theory of permutation group in my paper, Ref. [24] or in some earlier papers.

From the analysis of the discussed above theorem follows that for different values of the total spin S, the expression for obtained electron density does not changed and have the same value as for wave function presented as a single Slater determinant. In this connection, it should be mentioned that about the ambiguity of the description by the electron density was known many years ago. I will cite two publications: more than 40 years ago Harriman [30] demonstrated that for each electron density $\rho$ it can be constructed an arbitrary number of orthonormal orbitals, while in 2001 Cappelle and Vignale [31] showed that at the LSDA approximation it can be constructed different sets of potentials having the same ground state density.

Mean-while, different methods allowing taking into account the spin multiplet structure were developed, see Ref. [24] and references therein. In next section, I will discuss two groups of these methods.

3. Methods allowing take into account the spin multiplet structure within the DFT approach

It should be mentioned that the methods, taking into account the spin multiplet structure, are beyond the KS formalism. In most of these methods the $\mathrm{\Psi }$-formalism is used. I will consider two widely used approaches:

1. One of the first publications, in which the spin multiplet structure was taken into account in the frame of DFT, was the paper by Ziegler et al. [32]. For each value of the total spin S, they built the appropriate combinations of the Slater determinants. The factor giving the value of spin was obtained by correcting the exchange energy, $E_X$. In publications [33, 34, 35], in which scientists applied the Ziegler et al. approach, only the exchange energy was considered. Mineva et al. [36] are stressed that the scientists [33, 34, 35] developing after Ziegler et al. the methods that allows to take into account the conception of spin in the DFT studies, are considered only the exchange energy, $E_X$, and did not consider the correlation energy, $E_C$. This led to the incorrect multiplet structure. To the best of my knowledge, this drawback has not been discussed in the DFT community.

   In some applications, the Ziegler et al. method was named as Multiplet Structure Method, or shortly MSM. In the following text, I will use this abbreviation. Usually the exchange-correlation functional, $E_{\mathit{XC}}$, is presented as a sum of exchange, $E_X$, and correlation $E_C$, energies:

   $$E_{\mathit{xc}}=E_x+E_c.$$E9

   In the studies based or developing MSM approach [32, 33, 34, 35], the value of the total spin S was found using only the exchange energy $E_X$, and then applied to the total $E_{\mathit{xc}}$. It is evident that the contribution of the correlation energy $E_X$, is not the same as the exchange energy. Therefore, if $E_c$ is multiplied by the same factor as it was found for $E_x$ then it will give a wrong multiplet structure, because the exchange and correlations functionals should have a different dependence on S.

2. In the second group of methods [37, 38, 39, 40, 41], designated as restricted open-shell Kohn-Sham (ROKS) method, the open-shell theory of Roothaan [42] was used. In the first publication by Russo et al. [37], the Hamiltonian of Roothaan [42] was used, but the exchange term was replaced by the exchange-correlation functional. The authors [38, 39, 40, 41], combined the ROKS methods with the MSM approach. As a result, the methods elaborated in publications [38, 39, 40, 41] carried the same mistakes as the first ROKS method created by Russo et al. [37]. They do not provide the correct spin value for the correlation functional.

As stated above, both approaches MSM and ROKS, do not provide the correct total spin S for the correlation functionals.

According to calculations by Illas et al. [43], the ferromagnetic coupling is exaggerated, if the DFT method approaches are used. This agreed with our precise Mn₂ calculations performed with Mavridis group [44]. In most of calculations, Mn₂ had the ferromagnetic ground state with maximum value of the total spin, $S=5$. I would like also to mention that in the following article [45] Illas and collaborators, using the Filatov-Shaik ROKS method [39, 40], failed to improve the agreement with experiment. The reason is that, as we noted above, in the ROKS calculation, the correlation energy was not considered.

The definition of the correlation energy was given by Löwdin [46] many years ago. According to it:

$$E_{\mathit{\text{corr}}}=E_{\mathit{\text{exact}}}+E_{\mathit{HF}}.$$E10

The exact quantum-mechanical calculations can be performed only for small electron systems, for larger systems, the “exact” energy will depend upon the method used for its calculation. Thus, the correlation energy is method dependent. It should be also mentioned that the correlation energy has not an analytical expression, which leads to some problems in its applications, see Ref. [24].

As it was noted in Introduction, the total spin S of the system can be introduced only at the two-particle reduced density matrix level. The modern state of the development of the two-particle reduced density matrix formalism was discussed in a large number of articles [47, 48, 49, 50, 51, 52, 53, 54, 55]. Unfortunately, the spin problem still has not been considered by the DFT community.

4. Symmetry properties of the density matrix; degenerate states

In an elegant proof, Hohenberg and Kohn [1] laid down the theoretical foundation of the DFT theory. In their fundamental paper, the degeneracy was not treated, since they considered the ground state, which very rare is degenerated, as it is in the case of O₂ molecule.

Very soon, in the DFT community it was accepted that the Levy-Lieb [56, 57] constraint search procedure allows to study the degenerate states in the DFT calculations. First, it was shown by Levy [56] in 1969 and then in 1983 by Lieb [57], who applied more abstract mathematical approach. I would like to mention that Bersuker [58] was the first who criticized the possibility of application of the DFT approach to degenerate states. Bersuker considered it on the special case of the Jahn-Teller effect. Let us mention that according to the following from quantum mechanics the Born-Oppenheimer approximation (the molecules can be calculated only at this approximation) the vibronic interaction mixed the electronic and nuclear motions, and the electron and nuclear densities may not be constructed. Thus, the Levy-Lieb [56, 57] constraint search procedure contradicts quantum mechanics.

In Section 2, we already discussed that the author proved, see Ref. [24], the theorem that the electron density of the arbitrary N-electron system, defined in Eq. (1), does not depend upon the total spin S and always preserves the same form as it is for a single-determinantal wave function. From this theorem follows that the wave function of N-electron system does not depend on the degeneracy of the state and on its symmetry as well. It was proved using the permutation group apparatus, described in detail in chapter 2 of my book [28], see also Appendix to this chapter.

In general, it can be two types of degenerate states, the spatial and the spin degeneracy. In the case of the spatial degeneracy, the particles are described by spatial wave functions, although when they are degenerated in the spin space, they are described by spin wave functions. For constructing the degenerate in space wave function, the point group formalism should be used. Let us consider a point group G with $g$ elements. In the book [28], the author constructed the wave functions belonging to the $f_{\alpha }$-dimensional representation ${Г}^{\left(\alpha \right)}$ of an arbitrary point group $\mathit{G}$ as:

$${\mathrm{\Psi }}_{\mathit{ik}}^{\left(\alpha \right)}=\frac{f_{\alpha }}{g}\sum _R{\mathrm{\Gamma }}_{\mathit{ik}}^{\left(\alpha \right)}{\left(R\right)}^{\ast }R{\mathrm{\Psi }}_0,$$E11

where ${\mathrm{\Gamma }}_{\mathit{ik}}^{\left(\alpha \right)}\left(R\right)$ are the matrix elements of the representation ${Г}^{\left(\alpha \right)}$ and the sum in Eq. (11) is taken over all $g$ elements of the group G. The operations R of the group G are acting on some non-symmetrized product ${\mathrm{\Psi }}_0$ of one-electron orbitals. If in ${\mathrm{\Psi }}_{\mathit{ik}}^{\left(\alpha \right)}$ the second index $k$ is fixed, then $f_{\alpha }$ function form a basis for the representation ${Г}^{\left(\alpha \right)}$, each different indices $k$ enumerates different bases.

If we have a $f_{\alpha }$ degenerate state, each of its $f_{\alpha }$ functions describe the system with the same probability and a pure state (the state described by wave function) cannot be selected. Therefore, the degenerate state must be considered as a mixed state, in which its basis functions enter the density with the same weight factors. The diagonal element of the density matrix in the case of degenerate state is written as:

$$D_k^{\left(\alpha \right)}=\frac{1}{f_{\alpha }}\sum _{i=1}^{f_{\alpha }}{\left|{\mathrm{\Psi }}_{\mathit{ik}}^{\left(\alpha \right)}\right|}^2.$$E12

Using expression (12), I proved the theorem [24], according to which, the diagonal elements of the full density matrix is invariant for all operations of the group symmetry of the state, that is, it is a group invariant. It was proved that for every operation $R$ of group $G$ and all its irreducible representations ${Г}^{\left(\alpha \right)}$

$${\mathit{RD}}_k^{\alpha }=D_k^{\alpha }.$$E13

This means that the diagonal element of the full density matrix (and all reduced density matrices as well) transforms according to the totally symmetric one-dimensional representation $A_1$ of $G$ regardless of the dimension of representation ${Г}^{\left(\alpha \right)}$. It was proved for an arbitrary point group, but it is correct for any finite group. For the permutation group, this result was used in my publications [24, 59, 60] in analysis of the foundations of PEP. In these articles, I analyzed the case when PEP is not fulfilled and except of symmetrical and antisymmetrical states, an arbitrary permutation symmetry, including degenerate permutation states, are permitted. I have showed that if PEP is not fulfilled, this leads to contradictions with the concepts of particle independence and their identity. It was rigorously proved that the particles, described by wave functions with the permutation symmetry not allowed by PEP, may not exist in our Nature.

The arguments presented in Refs. [24, 59, 60], see also book [5], can be considered as a theoretical substantiation of PEP. They explained why in our Nature only completely symmetric or antisymmetric states, corresponding to one-dimensional representations of the permutation group, are realized. From this result, the important consequence follows

We may not expect that in future some unknown elementary particles can be discovered that are not fermions or bosons.

On the other hand, according to the so-called fractional statistics, which is valid in the 2D-space, a continuum of intermedium cases between boson and fermion particles can exist, see subsection 5.4 in book [5]. As was showed by Leinaas and Myrheim [61] in their pioneer paper, in 2D-space can exist a continuum of states between boson and fermion symmetry. After Leinaas and Myrheim [61], Wilczek [62] introduced in 2D-space the anyons, which obey any statistics. However, we should take into account that anyons are quasiparticles defined in 2D-space. The real particles can exist only in 3D-space, and according to PEP, formulated for all elementary particles, see Eq. (7), the elementary particles can obey only the boson or fermion symmetry. It is important to stress that the discovery of the fractional statistics does not contradict PEP.

All experimental data, see my recent review [63], confirm the Pauli Exclusion Principle. Different very precise experiments did not show any Pauli-forbidden transitions.

This is confirmed also by very precise calculations of H₂ molecule [7], in which, certainly, PEP was taken into account. The quantum mechanical calculations of the H₂ dissociation energy and its first ionization potential [64, 65] are in a complete agreement with very precise experimental values, see Table 1.1 in [7]. From this follows not only an additional confirmation of PEP, but also a rather general conclusion that molecules obey the same quantum-mechanical laws that obey traditionally physical objects: atoms and solids; at nanoscale we should not distinguish between chemical and physical systems.

In the end, I would like to note that in some papers the authors claimed that they developed the non-Born-Oppenheimer DFT in the frame of the electron density approach. These publications were analyzed in my first paper on DFT limits [24], where it was shown that in spite of the authors claims, their formalisms must be attributed to the Born-Oppenheimer approximation.

Acknowledgments

I am grateful to Maestro Ronald Columbié-Leyva for new references and for the technical and software support.

Conflict of interest

No potential conflict of interest was reported by the author.

The permutation symmetry is classified according to the irreducible representations of the permutation group ${\mathit{\pi }}_{\mathit{N}}$.¹ The latter are labeled by the Young diagrams

$$\left[\lambda \right]=\left[{\lambda }_1{\lambda }_2\dots {\lambda }_k\right],$$

$${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_k,\sum _{i=1}^k{\lambda }_i=N.$$E14

where ${\lambda }_i$ is represented by a row of ${\lambda }_i$ cells. The presence of several rows of equal length ${\lambda }_i$ is convenient to indicate by a power of ${\lambda }_i$. For example,

At present, the apparatus of permutations groups cannot be described without using the Young diagrams, I would like to note here some unusual details of the biography of Alfred Young. He was a country clergyman and has not any mathematical education. Young published studies were extending from 1900 to 1935, and in total he published 8 papers. The keystone of his studies was the reduction of the permutation groups to its irreducible representations in an explicit form. It is quite remarkable the gap of 25 years between his second paper in 1902 and the third in 1927. This gap will not be surprising, if we take into account that Young was a clergyman with numerous clerical duties.

It is obvious that one can form from two cells only two Young diagrams:

For the permutation group of three elements, ${\mathit{\pi }}_{\mathit{3}}$, one can form from three cells three Young diagrams:

The group ${\mathit{\pi }}_{\mathit{4}}$ has five Young diagrams:

As we mentioned in the beginning of this Appendix, the representation of the permutation group ${\mathit{\pi }}_{\mathit{N}}$ are labeled by the Young diagram $\left[\lambda \right]$. The rules how to find from a given Young diagram the matrices of representation of the permutation group were formulated by Young, but they were very complex. The Japanese mathematician Yamanouchi considerably simplified these rules. At present, it is called Young-Yamanouchi representation, which is described in detail in my book [28].

For each irreducible representation ${\mathrm{\Gamma }}^{\left[\lambda \right]}$, the normalized basis functions can be easily constructed by the Young operator,

$${\omega }_{\mathit{rt}}^{\left[\lambda \right]}=\sqrt{\frac{f_{\lambda }}{N!}}\sum _P{\mathrm{\Gamma }}_{\mathit{rt}}^{\left[\lambda \right]}\left(P\right)P.$$E15

In Eq. (15), ${\mathrm{\Gamma }}_{\mathit{rt}}^{\left[\lambda \right]}\left(P\right)$ are the matrix elements of the representation ${\mathrm{\Gamma }}^{\left[\lambda \right]}$, $f_{\lambda }$ is the dimension of ${\mathrm{\Gamma }}^{\left[\lambda \right]}$ and the sum over operations $P$ are taken for all $N!$ permutations of the group π_N. The normalized functions are obtained by acting the operator (15) on some non-symmetrized product of one-electron orbitals,

$${\mathrm{\Phi }}_0={\phi }_1\left(1\right){\phi }_2\left(2\right)\dots {\phi }_N\left(N\right).$$E16

The normalized functions

$${\mathrm{\Phi }}_{\mathit{rt}}^{\left[\lambda \right]}={\omega }_{\mathit{rt}}^{\left[\lambda \right]}{\mathrm{\Phi }}_0=\sqrt{\frac{f_{\lambda }}{N!}}\sum _P{\mathrm{\Gamma }}_{\mathit{rt}}^{\left[\lambda \right]}\left(P\right)P{\mathrm{\Phi }}_0,$$E17

are transformed in accordance with the representation ${\mathrm{\Gamma }}^{\left[\lambda \right]}$. If in ${\mathrm{\Phi }}_{\mathit{rt}}^{\left[\lambda \right]}$, Eq. (17), the second index $t$ is fixed, then $f_{\lambda }$ function form a basis for the representation ${\mathrm{\Gamma }}^{\left[\lambda \right]}$, each different index $k$ enumerates different bases. Let us prove this statement applying an arbitrary permutation $Q$ of the group ${\mathit{\pi }}_{\mathit{N}}$ to the function (17):

$$Q{\mathrm{\Phi }}_{\mathit{rt}}^{\left[\lambda \right]}=\sqrt{\frac{f_{\lambda }}{N!}}\sum _P{\mathrm{\Gamma }}_{\mathit{rt}}^{\left[\lambda \right]}\left(P\right)\mathit{QP}{\mathrm{\Phi }}_0=\sqrt{\frac{f_{\lambda }}{N!}}\sum _P{\mathrm{\Gamma }}_{\mathit{rt}}^{\left[\lambda \right]}\left(P\right)\left(Q^{-1}R\right)R{\mathrm{\Phi }}_0.$$E18

Using the invariance properties of a sum over all group elements and the property of orthogonal matrices, we obtain the matrix element of the product of permutations as a product of matrix elements.

$${\mathrm{\Gamma }}_{\mathit{rt}}^{\left[\lambda \right]}\left(P\right)\left(Q^{-1}R\right)=\sum _u{\mathrm{\Gamma }}_{\mathit{ru}}^{\left[\lambda \right]}\left(Q^{-1}\right){\mathrm{\Gamma }}_{\mathit{ut}}^{\left[\lambda \right]}\left(R\right)=\sum _u{\mathrm{\Gamma }}_{\mathit{ur}}^{\left[\lambda \right]}\left(Q\right){\mathrm{\Gamma }}_{\mathit{ut}}^{\left[\lambda \right]}\left(R\right).$$E19

In Eq. (19) we denoted the permutation $\mathit{QP}$ by $R$. Substituting (19) in (18), we obtain finally

$$Q{\mathrm{\Phi }}_{\mathit{rt}}^{\left[\lambda \right]}=\sqrt{\frac{f_{\lambda }}{N!}}\sum _u{\mathrm{\Gamma }}_{\mathit{ur}}^{\left[\lambda \right]}\left(Q\right)\left(\sum _R{\mathrm{\Gamma }}_{\mathit{ut}}^{\left[\lambda \right]}\left(R\right)R{\mathrm{\Phi }}_0\right)=\sum _u{\mathrm{\Gamma }}_{\mathit{ur}}^{\left[\lambda \right]}\left(Q\right){\mathrm{\Phi }}_{\mathit{ut}}^{\left[\lambda \right]}.$$E20

Thus, if the second index $t$ is fixed, then $f_{\lambda }$ functions form a basis for the representation ${\mathrm{\Gamma }}^{\left[\lambda \right]}$, and each different index $k$ enumerates different bases. It is also important to mention that in the function ${\mathrm{\Phi }}_{\mathit{rt}}^{\left[\lambda \right]}$, index $r$ characterizes the symmetry under permutation of the arguments, while index $t$ characterizes the symmetry under permutation of the one-particle functions ${\phi }_{\mathrm{a}}$.