Pair probability density
Under Hartree-Fock , the wave function is written as Slater determinant:
Ψ A S = ∣ χ 1 ( x 1 ) χ 2 ( x 1 ) ⋯ χ N ( x 1 ) χ 1 ( x 2 ) χ 2 ( x 2 ) ⋯ χ N ( x 2 ) ⋮ ⋮ ⋮ χ 1 ( x N ) χ 2 ( x N ) ⋯ χ N ( x N ) ∣ , (1) \Psi_{A S}=\left|\begin{array}{cccc}
\chi_{1}\left(\mathbf{x}_{1}\right) & \chi_{2}\left(\mathbf{x}_{1}\right) & \cdots & \chi_{N}\left(\mathbf{x}_{1}\right) \\
\chi_{1}\left(\mathbf{x}_{2}\right) & \chi_{2}\left(\mathbf{x}_{2}\right) & \cdots & \chi_{N}\left(\mathbf{x}_{2}\right) \\
\vdots & \vdots & & \vdots \\
\chi_{1}\left(\mathbf{x}_{N}\right) & \chi_{2}\left(\mathbf{x}_{N}\right) & \cdots & \chi_{N}\left(\mathbf{x}_{N}\right)
\end{array}\right|,
\tag{1} Ψ A S = χ 1 ( x 1 ) χ 1 ( x 2 ) ⋮ χ 1 ( x N ) χ 2 ( x 1 ) χ 2 ( x 2 ) ⋮ χ 2 ( x N ) ⋯ ⋯ ⋯ χ N ( x 1 ) χ N ( x 2 ) ⋮ χ N ( x N ) , ( 1 ) note that the normalization factor $(N!)^{-1/2}$ is ignored here so that $\braket{\Psi_{A S}|\Psi_{A S}}=N!$ where $N$ is the number of electron in the system.
Let’s first consider two electron case. If the two electrons have the same spin , then we can explicitly write out the Slater determinant as:
Ψ A S = ∣ χ 1 ( x 1 ) χ 2 ( x 1 ) χ 1 ( x 2 ) χ 2 ( x 2 ) ∣ , \Psi_{A S}=\left|\begin{array}{cc}
\chi_{1}\left(\mathbf{x}_{1}\right) & \chi_{2}\left(\mathbf{x}_{1}\right) \\
\chi_{1}\left(\mathbf{x}_{2}\right) & \chi_{2}\left(\mathbf{x}_{2}\right) \\
\end{array}\right|, Ψ A S = χ 1 ( x 1 ) χ 1 ( x 2 ) χ 2 ( x 1 ) χ 2 ( x 2 ) , where $\chi_1 = \psi_1 \otimes \beta_1$. $\psi$ signifies the spacial orbitals and $\beta$ is is the spin orbital of the down spin channel (or for the other spin channel $\alpha$).
Here, all $\chi$ are with $\beta$ spin for we are going to consider two same spin electrons.
For spin orbitals we have the orthornormal relations:
∫ α ( w ) β ( w ) d w = 0 ∫ β ( w ) β ( w ) d w = 1 \begin{aligned}
\int \alpha(w) \beta(w) dw &= 0 \\
\int \beta(w) \beta(w) dw &= 1
\end{aligned} ∫ α ( w ) β ( w ) d w ∫ β ( w ) β ( w ) d w = 0 = 1 If we expand out the determinant:
Ψ A S = [ ψ 1 ( r 1 ) β ( ω 1 ) ψ 2 ( r 2 ) β ( ω 2 ) − ψ 1 ( r 2 ) β ( ω 2 ) ψ 2 ( r 1 ) β ( ω 1 ) ] . \Psi_{A S}= [\psi_1(r_1)\beta(\omega_1)\psi_2(r_2)\beta(\omega_2)-\psi_1(r_2)\beta(\omega_2)\psi_2(r_1)\beta(\omega_1)]. Ψ A S = [ ψ 1 ( r 1 ) β ( ω 1 ) ψ 2 ( r 2 ) β ( ω 2 ) − ψ 1 ( r 2 ) β ( ω 2 ) ψ 2 ( r 1 ) β ( ω 1 )] . The probability of finding the first electron at $r_1$ with $dr_1$ and the second electron at $r_2$ with $dr_2$ is:
∣ Ψ ∣ 2 d r 1 d r 2 d ω 1 d ω 2 = ∣ ( ψ 1 ( r 1 ) β ( ω 1 ) ψ 2 ( r 2 ) β ( ω 2 ) − ψ 1 ( r 2 ) β ( ω 2 ) ψ 2 ( r 1 ) β ( ω 1 ) ) ∣ 2 d r 1 d r 2 d ω 1 d ω 2 |\Psi|^2 dr_1 dr_2 d\omega_1 d\omega_2 = |(\psi_1(r_1)\beta(\omega_1)\psi_2(r_2)\beta(\omega_2)-\psi_1(r_2)\beta(\omega_2)\psi_2(r_1)\beta(\omega_1))|^2 dr_1 dr_2 d\omega_1 d\omega_2 ∣Ψ ∣ 2 d r 1 d r 2 d ω 1 d ω 2 = ∣ ( ψ 1 ( r 1 ) β ( ω 1 ) ψ 2 ( r 2 ) β ( ω 2 ) − ψ 1 ( r 2 ) β ( ω 2 ) ψ 2 ( r 1 ) β ( ω 1 )) ∣ 2 d r 1 d r 2 d ω 1 d ω 2 Integrating out $\omega_1$ and $\omega_2$, we get:
P ( r 1 , r 2 ) d r 1 d r 2 = [ ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 + ∣ ψ 1 ( r 2 ) ∣ 2 ∣ ψ 2 ( r 1 ) ∣ 2 ] − ψ 1 ∗ ( r 1 ) ψ 2 ( r 1 ) ψ 2 ∗ ( r 2 ) ψ 1 ( r 2 ) − ψ 1 ( r 1 ) ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ψ 1 ∗ ( r 2 ) d r 1 d r 2 , (2) \begin{aligned}
P(r_1,r_2) dr_1 dr_2 = & [|\psi_1(r_1)|^2 |\psi_2(r_2)|^2 + |\psi_1(r_2)|^2 |\psi_2(r_1)|^2]\\ & - \psi^*_1(r_1)\psi_2(r_1)\psi^*_2(r_2)\psi_1(r_2)\\
& - \psi_1(r_1)\psi^*_2(r_1)\psi_2(r_2)\psi^*_1(r_2) dr_1 dr_2, \tag{2}
\end{aligned} P ( r 1 , r 2 ) d r 1 d r 2 = [ ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 + ∣ ψ 1 ( r 2 ) ∣ 2 ∣ ψ 2 ( r 1 ) ∣ 2 ] − ψ 1 ∗ ( r 1 ) ψ 2 ( r 1 ) ψ 2 ∗ ( r 2 ) ψ 1 ( r 2 ) − ψ 1 ( r 1 ) ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ψ 1 ∗ ( r 2 ) d r 1 d r 2 , ( 2 ) note that we’ve taken advantage of the orthornormal relation bettween $\alpha$ and $\beta$.
Eq. 2 gives us the probability of finding two same spin electrons at $r_1$ and $r_2$ (within the volume of $dr_1$ and $dr_2$), we call this probability the pair probability.
It is sometimes also called the second-order density matrix (
Next, we can add change the look of this pair probability by adding some additional terms:
P ( r 1 , r 2 ) d r 1 d r 2 = [ ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 + ∣ ψ 1 ( r 2 ) ∣ 2 ∣ ψ 2 ( r 1 ) ∣ 2 − ψ 1 ∗ ( r 1 ) ψ 2 ( r 1 ) ψ 2 ∗ ( r 2 ) ψ 1 ( r 2 ) − ψ 1 ( r 1 ) ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ψ 1 ∗ ( r 2 ) ] d r 1 d r 2 = { ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 + ∣ ψ 1 ( r 2 ) ∣ 2 ∣ ψ 2 ( r 1 ) ∣ 2 + [ ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 1 ( r 2 ) ∣ 2 + ∣ ψ 2 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 − ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 1 ( r 2 ) ∣ 2 − ∣ ψ 2 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 ] − ψ 1 ∗ ( r 1 ) ψ 2 ( r 1 ) ψ 2 ∗ ( r 2 ) ψ 1 ( r 2 ) − ψ 1 ( r 1 ) ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ψ 1 ∗ ( r 2 ) } d r 1 d r 2 = { [ ∣ ψ 1 ( r 1 ) ∣ 2 + ∣ ψ 2 ( r 1 ) ∣ 2 ] ⋅ [ ∣ ψ 1 ( r 2 ) ∣ 2 + ∣ ψ 2 ( r 2 ) ∣ 2 ] − ∣ ψ 1 ∗ ( r 1 ) ψ 1 ( r 2 ) + ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ∣ 2 } d r 1 d r 2 = [ ∑ i ψ i ∗ ( r 1 ) ψ i ( r 1 ) ∑ j ψ j ∗ ( r 2 ) ψ j ( r 2 ) − ∣ ∑ k ψ k ∗ ( r 1 ) ψ k ( r 2 ) ∣ 2 ] d r 1 d r 2 = [ ρ ( r 1 ) ρ ( r 2 ) − ∣ ρ ( r 1 , r 2 ) ∣ 2 ] d r 1 d r 2 , (3) \begin{aligned}
P(r_1,r_2) dr_1 dr_2 &= [|\psi_1(r_1)|^2 |\psi_2(r_2)|^2 + |\psi_1(r_2)|^2 |\psi_2(r_1)|^2\\ & \quad - \psi^*_1(r_1)\psi_2(r_1)\psi^*_2(r_2)\psi_1(r_2) - \psi_1(r_1)\psi^*_2(r_1)\psi_2(r_2)\psi^*_1(r_2)] dr_1 dr_2 \\
&= \{ |\psi_1(r_1)|^2 |\psi_2(r_2)|^2 + |\psi_1(r_2)|^2 |\psi_2(r_1)|^2\\ & \textcolor{red}{\quad + [ |\psi_1(r_1)|^2 |\psi_1(r_2)|^2 + |\psi_2(r_1)|^2 |\psi_2(r_2)|^2} \textcolor{red}{ - |\psi_1(r_1)|^2 |\psi_1(r_2)|^2 - |\psi_2(r_1)|^2 |\psi_2(r_2)|^2 ]}\\
& \quad - \psi^*_1(r_1)\psi_2(r_1)\psi^*_2(r_2)\psi_1(r_2) - \psi_1(r_1)\psi^*_2(r_1)\psi_2(r_2)\psi^*_1(r_2) \} dr_1 dr_2\\
&= \{ [ |\psi_1(r_1)|^2+|\psi_2(r_1)|^2] \cdot [|\psi_1(r_2)|^2+|\psi_2(r_2)|^2]\\ & \quad-|\psi^*_1(r_1)\psi_1(r_2)+\psi^*_2(r_1)\psi_2(r_2)|^2 \} dr_1 dr_2\\
&= [ \sum_i \psi^*_i(r_1)\psi_i(r_1) \sum_j \psi^*_j(r_2)\psi_j(r_2) - |\sum_k \psi^*_k(r_1)\psi_k(r_2)|^2 ]dr_1 dr_2\\
&= [ \rho(r_1)\rho(r_2) -|\rho(r_1,r_2)|^2 ]dr_1 dr_2, \tag{3}
\end{aligned} P ( r 1 , r 2 ) d r 1 d r 2 = [ ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 + ∣ ψ 1 ( r 2 ) ∣ 2 ∣ ψ 2 ( r 1 ) ∣ 2 − ψ 1 ∗ ( r 1 ) ψ 2 ( r 1 ) ψ 2 ∗ ( r 2 ) ψ 1 ( r 2 ) − ψ 1 ( r 1 ) ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ψ 1 ∗ ( r 2 )] d r 1 d r 2 = { ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 + ∣ ψ 1 ( r 2 ) ∣ 2 ∣ ψ 2 ( r 1 ) ∣ 2 + [ ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 1 ( r 2 ) ∣ 2 + ∣ ψ 2 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 − ∣ ψ 1 ( r 1 ) ∣ 2 ∣ ψ 1 ( r 2 ) ∣ 2 − ∣ ψ 2 ( r 1 ) ∣ 2 ∣ ψ 2 ( r 2 ) ∣ 2 ] − ψ 1 ∗ ( r 1 ) ψ 2 ( r 1 ) ψ 2 ∗ ( r 2 ) ψ 1 ( r 2 ) − ψ 1 ( r 1 ) ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ψ 1 ∗ ( r 2 )} d r 1 d r 2 = {[ ∣ ψ 1 ( r 1 ) ∣ 2 + ∣ ψ 2 ( r 1 ) ∣ 2 ] ⋅ [ ∣ ψ 1 ( r 2 ) ∣ 2 + ∣ ψ 2 ( r 2 ) ∣ 2 ] − ∣ ψ 1 ∗ ( r 1 ) ψ 1 ( r 2 ) + ψ 2 ∗ ( r 1 ) ψ 2 ( r 2 ) ∣ 2 } d r 1 d r 2 = [ i ∑ ψ i ∗ ( r 1 ) ψ i ( r 1 ) j ∑ ψ j ∗ ( r 2 ) ψ j ( r 2 ) − ∣ k ∑ ψ k ∗ ( r 1 ) ψ k ( r 2 ) ∣ 2 ] d r 1 d r 2 = [ ρ ( r 1 ) ρ ( r 2 ) − ∣ ρ ( r 1 , r 2 ) ∣ 2 ] d r 1 d r 2 , ( 3 ) where we have used the one-body density matrix:
ρ ( r 1 ) = ρ ( r 1 , r 1 ) = ∑ i ψ i ∗ ( r 1 ) ψ i ( r 1 ) ρ ( r 1 , r 2 ) = ∑ i ψ i ∗ ( r 1 ) ψ i ( r 2 ) (4) \begin{aligned}
\rho(r_1) &= \rho(r_1, r_1) = \sum_i \psi^*_i(r_1)\psi_i(r_1) \\
\rho(r_1, r_2) &= \sum_i \psi^*_i(r_1)\psi_i(r_2) \tag{4}
\end{aligned} ρ ( r 1 ) ρ ( r 1 , r 2 ) = ρ ( r 1 , r 1 ) = i ∑ ψ i ∗ ( r 1 ) ψ i ( r 1 ) = i ∑ ψ i ∗ ( r 1 ) ψ i ( r 2 ) ( 4 ) The one body density matrix also have the following property:
∫ ∣ ρ 1 ( r 1 , r 2 ) ∣ 2 d r 2 = ρ ( r 1 ) . (6) \int\left|\rho_{1}(r_{1},r_{2})\right|^{2} dr_{2}=\rho(r_{1}). \tag{6} ∫ ∣ ρ 1 ( r 1 , r 2 ) ∣ 2 d r 2 = ρ ( r 1 ) . ( 6 ) Electron Localization Function (ELF)
If an electron (with spin $\sigma$) is located with certainty at position $r_{1}$, then the conditional probability of finding a second like-spin electron at position $r_{2}$ is obtained by dividing the pair probability by the probability of finding the first electron at $r_{1}$, i.e. the total charge density at $r_{1}$ (see the definition on condictional probability ):
P cond σ σ ( r 1 , r 2 ) = P σ σ ( r 1 , r 2 ) ρ σ ( r 1 ) = ρ σ ( r 2 ) − ∣ ρ σ ( r 1 , r 2 ) ∣ 2 ρ σ ( r 1 ) , (7) \begin{aligned}
P^{\sigma\sigma}_\text{cond}\left(r_{1}, r_{2}\right) &= \frac{ P^{\sigma\sigma}\left(r_{1}, r_{2}\right)}{\rho^{\sigma}(r_{1})}\\
&= \rho^{\sigma}(r_{2}) - \frac{\left|\rho^{\sigma}(r_{1},r_{2})\right|^2}{\rho^{\sigma}(r_{1})}, \tag{7}
\end{aligned} P cond σσ ( r 1 , r 2 ) = ρ σ ( r 1 ) P σσ ( r 1 , r 2 ) = ρ σ ( r 2 ) − ρ σ ( r 1 ) ∣ ρ σ ( r 1 , r 2 ) ∣ 2 , ( 7 ) using the properties of $\rho^{\sigma \sigma}(r_{1},r_{2})$, we can see that,
P cond σ σ ( r 1 , r 1 ) = 0 , (8) P^{\sigma\sigma}_{\text {cond }}(r_{1},r_{1})=0, \tag{8} P cond σσ ( r 1 , r 1 ) = 0 , ( 8 ) and (see Eq. 10 in becke1983 ),
∫ ∣ ρ σ ( r 1 , r 2 ) ∣ 2 ρ σ ( r 1 ) d r 2 = 1 , \int \frac{\left|\rho^{\sigma}(r_{1},r_{2})\right|^2}{\rho^{\sigma}(r_{1})} dr_2 = 1, ∫ ρ σ ( r 1 ) ∣ ρ σ ( r 1 , r 2 ) ∣ 2 d r 2 = 1 , so that,
∫ P cond σ σ ( r 1 , r 2 ) d r 2 = N − 1 , (9) \int P^{\sigma\sigma}_{\text {cond }}(r_{1},r_{2}) d r_{2}=N-1, \tag{9} ∫ P cond σσ ( r 1 , r 2 ) d r 2 = N − 1 , ( 9 ) where where $N$ is the total number electrons (with the same spin) in the system.
i.e., if a $\sigma$ spin electron is definitely at point $r_1$, then the total probability of finding another $\sigma$-spin electron elsewhere in the system is $N - 1$ (we have $N-1$ electrons left with the same spin, remember in the beginning we set spin of all electrons to be the same $\beta$). In our case, $N=2$ so the total probability of finding another $\sigma$-spin electron is $1$.
Now, if we take the the second term from Eq. 7, we see that it has the unit of a charge density.
This, turns out to be the so-called the exchange charge density (See becke1983 ):
ρ x σ σ ( r 1 , r 2 ) = ∣ ρ σ ( r 1 , r 2 ) ∣ 2 ρ σ ( r 1 ) . (10) \rho_{x}^{\sigma \sigma} (r_1,r_2) = \frac{\left|\rho^{\sigma}(r_{1},r_{2})\right|^2}{\rho^{\sigma}(r_{1})}. \tag{10} ρ x σσ ( r 1 , r 2 ) = ρ σ ( r 1 ) ∣ ρ σ ( r 1 , r 2 ) ∣ 2 . ( 10 ) and the numerator is called the exchange density.
Remember the Taylor expansion of $e^x$ around $x=0$ has the form of:
e x = ∑ n = 0 ∞ x n n ! = 1 + x + x 2 2 ! + x 3 3 ! + ⋯ , e^{x}=\sum_{n=0}^{\infty} \frac{x^{n}}{n !}=1+x+\frac{x^{2}}{2 !}+\frac{x^{3}}{3 !}+\cdots, e x = n = 0 ∑ ∞ n ! x n = 1 + x + 2 ! x 2 + 3 ! x 3 + ⋯ , comparing this to the definition of Taylor series:
f ( x ) = f ( a ) + f ′ ( a ) 1 ! ( x − a ) + f ′ ′ ( a ) 2 ! ( x − a ) 2 + f ′ ′ ′ ( a ) 3 ! ( x − a ) 3 + ⋯ , f(x) = f(a)+\frac{f^{\prime}(a)}{1 !}(x-a)+\frac{f^{\prime \prime}(a)}{2 !}(x-a)^{2}+\frac{f^{\prime \prime \prime}(a)}{3 !}(x-a)^{3}+\cdots, f ( x ) = f ( a ) + 1 ! f ′ ( a ) ( x − a ) + 2 ! f ′′ ( a ) ( x − a ) 2 + 3 ! f ′′′ ( a ) ( x − a ) 3 + ⋯ , or, if we switch to another coordinate system that’s centered around $a$:
f ( a + s ) = f ( a ) + f ′ ( a ) 1 ! ( s ) + f ′ ′ ( a ) 2 ! ( s ) 2 + f ′ ′ ′ ( a ) 3 ! ( s ) 3 + ⋯ , f(a+s) = f(a)+\frac{f^{\prime}(a)}{1 !}(s)+\frac{f^{\prime \prime}(a)}{2 !}(s)^{2}+\frac{f^{\prime \prime \prime}(a)}{3 !}(s)^{3}+\cdots, f ( a + s ) = f ( a ) + 1 ! f ′ ( a ) ( s ) + 2 ! f ′′ ( a ) ( s ) 2 + 3 ! f ′′′ ( a ) ( s ) 3 + ⋯ , we see that we can express the function at a specific point ($a+s$) as:
f ( a + s ) = e s d d x f ( x ) ∣ x = a = { 1 + s d d x + [ s d d x ] 2 2 ! + [ s d d x ] 3 3 ! + ⋯ } f ( a ) = f ( a ) + f ′ ( a ) 1 ! s + f ′ ′ ( a ) 2 ! s 2 + f ′ ′ ′ ( a ) 3 ! s 3 + ⋯ . (11) \begin{aligned}
f(a+s) &= e^{s\frac{d}{d x}} f(x)|_{x=a}\\
&= \left \{1 + s\frac{d}{d x} +\frac{[s\frac{d}{d x}]^{2}}{2 !}+\frac{[s\frac{d}{d x}]^{3}}{3 !}+\cdots \right \} f(a) \\
&= f(a)+\frac{f^{\prime}(a)}{1 !}s+\frac{f^{\prime \prime}(a)}{2 !}s^{2}+\frac{f^{\prime \prime \prime}(a)}{3 !}s^{3}+\cdots. \tag{11}
\end{aligned} f ( a + s ) = e s d x d f ( x ) ∣ x = a = { 1 + s d x d + 2 ! [ s d x d ] 2 + 3 ! [ s d x d ] 3 + ⋯ } f ( a ) = f ( a ) + 1 ! f ′ ( a ) s + 2 ! f ′′ ( a ) s 2 + 3 ! f ′′′ ( a ) s 3 + ⋯ . ( 11 ) Using Eq. 11, we can then Taylor expand the exchange charge density $\rho_{x}^{\sigma \sigma} (r_1,r_2)$ for $r_2$ near $r_1$:
ρ x σ σ ( r ⃗ 1 , r ⃗ 2 ) = ρ x σ σ ( r ⃗ 1 , r ⃗ 2 − r ⃗ 1 + r ⃗ 1 ) = ρ x σ σ ( r ⃗ 1 , s ⃗ + r ⃗ 1 ) ρ x σ σ ( r ⃗ 1 , r ⃗ 1 + s ⃗ ) = e s ⃗ ⋅ ∇ ⃗ r ⃗ ρ x σ σ ( r ⃗ 1 , r ⃗ ) ∣ r ⃗ → r ⃗ 1 , (12) \begin{aligned}
\rho_{x}^{\sigma \sigma} (\vec r_1,\vec r_2) &= \rho_{x}^{\sigma \sigma} (\vec r_1, \vec r_2- \vec r_1+ \vec r_1)\\
&= \rho_{x}^{\sigma \sigma} (\vec r_1,\vec s+\vec r_1)\\
\rho_{x}^{\sigma \sigma} (\vec r_1,\vec r_1 + \vec s) &= e^{\vec s \cdot \vec \nabla_{\vec r}} \rho_{x}^{\sigma \sigma} (\vec r_1,\vec r)|_{\vec r \to \vec r_1}, \tag{12}
\end{aligned} ρ x σσ ( r 1 , r 2 ) ρ x σσ ( r 1 , r 1 + s ) = ρ x σσ ( r 1 , r 2 − r 1 + r 1 ) = ρ x σσ ( r 1 , s + r 1 ) = e s ⋅ ∇ r ρ x σσ ( r 1 , r ) ∣ r → r 1 , ( 12 ) where we have defined $\vec s = \vec r_2 -\vec r_1$. As $\vec r_2 \to \vec r_1$, $\vec s \to 0$. Note that $r$ here is completely different from $r_2$ used before.
If we take the spherical average of the Taylor expansion $e^{\vec s \cdot \vec \nabla_{\vec r_2}}$ around $\vec r_2 \to \vec r_1$ with $\vec s$:
⟨ e s ⃗ ⋅ ∇ ⃗ r ⃗ 2 ⟩ = 1 4 π ∣ s ⃗ ∣ 2 ∫ ∫ e s ⃗ ⋅ ∇ ⃗ r ⃗ 2 ∣ s ⃗ ∣ 2 sin ( ϕ s ) d ϕ s d θ s \begin{aligned}
\langle e^{\vec s \cdot \vec \nabla_{\vec r_2}} \rangle &= \frac{1}{4\pi |\vec s|^2} \int \int e^{\vec s \cdot \vec \nabla_{\vec r_2}} |\vec s|^2 \sin(\phi_s)\, d\phi_s\, d \theta_s\\
\end{aligned} ⟨ e s ⋅ ∇ r 2 ⟩ = 4 π ∣ s ∣ 2 1 ∫∫ e s ⋅ ∇ r 2 ∣ s ∣ 2 sin ( ϕ s ) d ϕ s d θ s according to Eq. 12 in becke1983 , we arrived at:
⟨ e s ⃗ ⋅ ∇ ⃗ r ⃗ 2 ⟩ = sinh ( s ∇ r ⃗ 2 ) ( s ∇ r ⃗ 2 ) , (13) \langle e^{\vec s \cdot \vec \nabla_{\vec r_2}} \rangle =\frac{\sinh \left(s \nabla_{\vec r_2}\right)}{\left(s \nabla_{\vec r_2}\right)}, \tag{13} ⟨ e s ⋅ ∇ r 2 ⟩ = ( s ∇ r 2 ) sinh ( s ∇ r 2 ) , ( 13 ) Here, $\nabla_{\vec r_2}$ is the derivative along $\vec s = \vec r_2 - \vec r_1$ direction and is angle independent (it only acts on the coordinates of the second electron $\vec r_2$ and not $\vec s$ in the double integral, that’s why I haven’t wrote it as $\nabla_{\vec s}$, later on after I’ve applied this spherically averaged operator on to a state I’ll express it as $\nabla_{\vec s}$).
If we use the Taylor expansion of $\sinh$ function , we get:
⟨ e s ⃗ ⋅ ∇ ⃗ r ⃗ 2 ⟩ = 1 + 1 3 ! s 2 ∇ r ⃗ 2 2 + 1 5 ! s 4 ∇ r ⃗ 2 4 + 1 7 ! s 6 ∇ r ⃗ 2 6 + ⋯ , (14) \langle e^{\vec s \cdot \vec \nabla_{\vec r_2}} \rangle=1+\frac{1}{3 !} s^{2} \nabla_{\vec r_2}^{2}+\frac{1}{5 !} s^{4} \nabla_{\vec r_2}^{4}+\frac{1}{7 !} s^{6} \nabla_{\vec r_2}^{6}+\cdots, \tag{14} ⟨ e s ⋅ ∇ r 2 ⟩ = 1 + 3 ! 1 s 2 ∇ r 2 2 + 5 ! 1 s 4 ∇ r 2 4 + 7 ! 1 s 6 ∇ r 2 6 + ⋯ , ( 14 ) note that $s$ in Eq. 14 is a scalar.
Now we have the form of the spherically averaged Taylor expansion operator in Eq. 14, plugging this back into Eq. 12 gives us:
⟨ ρ x ( r ⃗ 1 , s ) ⟩ = ( 1 + 1 6 s 2 ∇ r ⃗ 2 2 + ⋯ ) ρ x ( r ⃗ 1 , r ⃗ 2 ) ∣ r ⃗ 2 → r ⃗ 1 = ( 1 + 1 6 s 2 ∇ s ⃗ 2 + ⋯ ) ρ x ( r ⃗ 1 , s ⃗ ) ∣ s ⃗ → 0 (15) \begin{aligned}
\left\langle\rho_{x}(\vec r_1, s)\right\rangle&=\left.\left(1+\frac{1}{6} s^{2} \nabla_{\vec r_2}^{2}+\cdots\right) \rho_{x}\left(\vec r_1, \vec r_2\right)\right|_{\vec r_2 \to \vec r_1}\\
&=\left.\left(1+\frac{1}{6} s^{2} \nabla_{\vec s}^{2}+\cdots\right) \rho_{x}\left(\vec r_1, \vec s\right)\right|_{\vec s \to 0}\\
\end{aligned} \tag{15} ⟨ ρ x ( r 1 , s ) ⟩ = ( 1 + 6 1 s 2 ∇ r 2 2 + ⋯ ) ρ x ( r 1 , r 2 ) r 2 → r 1 = ( 1 + 6 1 s 2 ∇ s 2 + ⋯ ) ρ x ( r 1 , s ) s → 0 ( 15 ) Cutting off the expansion up to second order, we need the result of operation with the Laplacian $\nabla_{\vec s}^2$ on the exchange charge $\rho_x$:
∇ r ⃗ 2 2 ρ x ( r ⃗ 1 , r ⃗ 2 ) ∣ r ⃗ 2 → r ⃗ 1 = ∇ s ⃗ 2 ∑ i ψ i ∗ ( r ⃗ 1 ) ψ i ( r ⃗ 2 ) ∑ j ψ j ( r ⃗ 1 ) ψ j ∗ ( r ⃗ 2 ) ∑ k ψ k ∗ ( r ⃗ 1 ) ψ k ( r ⃗ 1 ) ∇ s ⃗ 2 ρ x ( r ⃗ 1 , s ⃗ ) ∣ s ⃗ → 0 = ∇ s ⃗ 2 ∑ i ∑ j ψ i ∗ ( r ⃗ 1 ) ψ i ( r ⃗ 1 + s ⃗ ) ψ j ( r ⃗ 1 ) ψ j ∗ ( r ⃗ 1 + s ⃗ ) ∑ k ψ k ∗ ( r ⃗ 1 ) ψ k ( r ⃗ 1 ) = ∑ i ∑ j ψ i ∗ ( r ⃗ 1 ) ψ j ( r ⃗ 1 ) ∇ s ⃗ 2 [ ψ i ( r ⃗ 1 + s ⃗ ) ψ j ∗ ( r ⃗ 1 + s ⃗ ) ] ∑ k ψ k ∗ ( r ⃗ 1 ) ψ k ( r ⃗ 1 ) = ∑ i ∑ j [ ψ i ∗ ψ i ( ψ j ∇ 2 ψ j ∗ ) + ψ j ψ j ∗ ( ψ i ∗ ∇ 2 ψ i ) + ψ i ∗ ψ j ( 2 ∇ ψ i ∇ ψ j ∗ ) ] ∑ k ψ k ∗ ψ k . (16) \begin{aligned}
\nabla_{\vec r_2}^2 \rho_{x}\left(\vec r_1, \vec r_2\right)|_{\vec r_2 \to \vec r_1} &= \nabla_{\vec s}^2
\frac{\sum_i \psi^*_i(\vec r_1)\psi_i(\vec r_2) \sum_j \psi_j(\vec r_1)\psi^*_j(\vec r_2)}{\sum_k \psi^*_k(\vec r_1)\psi_k(\vec r_1)}\\
\nabla_{\vec s}^2 \rho_{x}\left(\vec r_1, \vec s\right)|_{\vec s \to 0} &= \nabla_{\vec s}^2\frac{\sum_i \sum_j \psi^*_i(\vec r_1)\psi_i(\vec r_1+\vec s) \psi_j(\vec r_1)\psi^*_j(\vec r_1+\vec s)}{\sum_k \psi^*_k(\vec r_1)\psi_k(\vec r_1)}\\
&= \frac{\sum_i \sum_j \psi^*_i(\vec r_1)\psi_j(\vec r_1) \nabla_{\vec s}^2 [\psi_i(\vec r_1+\vec s)\psi^*_j(\vec r_1+\vec s)]}{\sum_k \psi^*_k(\vec r_1)\psi_k(\vec r_1)}\\
&=\frac{\sum_i\sum_j \left [ \psi_i^* \psi_i (\psi_j \nabla^2 \psi_j^*) +\psi_j \psi_j^* (\psi_i^* \nabla^2 \psi_i)
+ \psi_i^* \psi_j (2 \nabla \psi_i \nabla \psi_j^*)
\right ]}
{\sum_k \psi^*_k\psi_k}.
\end{aligned} \tag{16} ∇ r 2 2 ρ x ( r 1 , r 2 ) ∣ r 2 → r 1 ∇ s 2 ρ x ( r 1 , s ) ∣ s → 0 = ∇ s 2 ∑ k ψ k ∗ ( r 1 ) ψ k ( r 1 ) ∑ i ψ i ∗ ( r 1 ) ψ i ( r 2 ) ∑ j ψ j ( r 1 ) ψ j ∗ ( r 2 ) = ∇ s 2 ∑ k ψ k ∗ ( r 1 ) ψ k ( r 1 ) ∑ i ∑ j ψ i ∗ ( r 1 ) ψ i ( r 1 + s ) ψ j ( r 1 ) ψ j ∗ ( r 1 + s ) = ∑ k ψ k ∗ ( r 1 ) ψ k ( r 1 ) ∑ i ∑ j ψ i ∗ ( r 1 ) ψ j ( r 1 ) ∇ s 2 [ ψ i ( r 1 + s ) ψ j ∗ ( r 1 + s )] = ∑ k ψ k ∗ ψ k ∑ i ∑ j [ ψ i ∗ ψ i ( ψ j ∇ 2 ψ j ∗ ) + ψ j ψ j ∗ ( ψ i ∗ ∇ 2 ψ i ) + ψ i ∗ ψ j ( 2∇ ψ i ∇ ψ j ∗ ) ] . ( 16 ) Note that we are Taylor expanding at $\vec r_2 \to \vec r_1$ ($\vec s \to 0$), thus, we have changed all $\vec r_2$ dependence to $\vec r_1$($\vec r_1+\vec s$).
Because∇ s ⃗ ψ i ( r ⃗ 1 + s ⃗ ) ∣ s ⃗ → 0 = ∇ r ⃗ 1 ψ i ( r ⃗ 1 ) \nabla_{\vec s} \psi_i(\vec r_1 + \vec s) |_{\vec s \to 0} = \nabla_{\vec r_1} \psi_i(\vec r_1) ∇ s ψ i ( r 1 + s ) ∣ s → 0 = ∇ r 1 ψ i ( r 1 )
The first two terms in Eq. 16 can be expressed as:
∑ i ∑ j [ ψ i ∗ ψ i ( ψ j ∇ 2 ψ j ∗ ) + ψ j ψ j ∗ ( ψ i ∗ ∇ 2 ψ i ) ] ∑ k ψ k ∗ ψ k = ∑ j ψ j ∇ 2 ψ j ∗ + ∑ i ψ i ∗ ∇ 2 ψ i = ∑ i ψ i ∇ 2 ψ i ∗ + ∑ i ψ i ∗ ∇ 2 ψ i + 2 ∑ i ∇ ψ i ∇ ψ i ∗ − 2 ∑ i ∇ ψ i ∇ ψ i ∗ = ∑ i ( ψ i ∇ 2 ψ i ∗ + 2 ∇ ψ i ∇ ψ i ∗ + ψ i ∗ ∇ 2 ψ i ) − 2 ∑ i ∇ ψ i ∇ ψ i ∗ , \begin{aligned}
&\frac{\sum_i\sum_j \left [ \psi_i^* \psi_i (\psi_j \nabla^2 \psi_j^*) +\psi_j \psi_j^* (\psi_i^* \nabla^2 \psi_i) \right ]}
{\sum_k \psi^*_k\psi_k}\\
&= \sum_j \psi_j \nabla^2 \psi_j^* + \sum_i \psi_i^* \nabla^2 \psi_i\\
&= \sum_i \psi_i \nabla^2 \psi_i^* + \sum_i \psi_i^* \nabla^2 \psi_i + 2 \sum_i \nabla \psi_i \nabla \psi_i^* - 2 \sum_i \nabla \psi_i \nabla \psi_i^*\\
&= \sum_i (\psi_i \nabla^2 \psi_i^* + 2 \nabla \psi_i \nabla \psi_i^* + \psi_i^* \nabla^2 \psi_i) - 2 \sum_i \nabla \psi_i \nabla \psi_i^*,
\end{aligned} ∑ k ψ k ∗ ψ k ∑ i ∑ j [ ψ i ∗ ψ i ( ψ j ∇ 2 ψ j ∗ ) + ψ j ψ j ∗ ( ψ i ∗ ∇ 2 ψ i ) ] = j ∑ ψ j ∇ 2 ψ j ∗ + i ∑ ψ i ∗ ∇ 2 ψ i = i ∑ ψ i ∇ 2 ψ i ∗ + i ∑ ψ i ∗ ∇ 2 ψ i + 2 i ∑ ∇ ψ i ∇ ψ i ∗ − 2 i ∑ ∇ ψ i ∇ ψ i ∗ = i ∑ ( ψ i ∇ 2 ψ i ∗ + 2∇ ψ i ∇ ψ i ∗ + ψ i ∗ ∇ 2 ψ i ) − 2 i ∑ ∇ ψ i ∇ ψ i ∗ , since
∇ ψ i ∗ = ( ∇ ψ i ) ∗ \nabla \psi_i^* = \left( \nabla \psi_i \right)^* ∇ ψ i ∗ = ( ∇ ψ i ) ∗ ∇ ψ i ∇ ψ i ∗ = ∣ ∇ ψ i ∣ 2 \nabla \psi_i \nabla \psi_{i}^{*} = \lvert \nabla \psi_i \rvert^{2} ∇ ψ i ∇ ψ i ∗ = ∣ ∇ ψ i ∣ 2
∑ i ∑ j [ ψ i ∗ ψ i ( ψ j ∇ 2 ψ j ∗ ) + ψ j ψ j ∗ ( ψ i ∗ ∇ 2 ψ i ) ] ∑ k ψ k ∗ ψ k = ∇ 2 ∑ i ∣ ψ i ∣ 2 − 2 ∑ i ∣ ∇ ψ i ∣ 2 = ∇ 2 ρ − 2 τ , (17) \begin{aligned}
&\frac{\sum_i\sum_j \left [ \psi_i^* \psi_i (\psi_j \nabla^2 \psi_j^*) +\psi_j \psi_j^* (\psi_i^* \nabla^2 \psi_i) \right ]}
{\sum_k \psi^*_k\psi_k}\\
&= \nabla^2 \sum_i |\psi_i|^2 - 2 \sum_i|\nabla \psi_i|^2\\
&= \nabla^2 \rho - 2\tau,
\end{aligned} \tag{17} ∑ k ψ k ∗ ψ k ∑ i ∑ j [ ψ i ∗ ψ i ( ψ j ∇ 2 ψ j ∗ ) + ψ j ψ j ∗ ( ψ i ∗ ∇ 2 ψ i ) ] = ∇ 2 i ∑ ∣ ψ i ∣ 2 − 2 i ∑ ∣∇ ψ i ∣ 2 = ∇ 2 ρ − 2 τ , ( 17 ) where we have used the positively defined kinetic energy density
τ = ∑ i ∣ ∇ ψ i ∣ 2 \tau=\sum_i|\nabla \psi_i|^2 τ = ∑ i ∣∇ ψ i ∣ 2
The third term in Eq. 16 is:
∑ i ∑ j [ ψ i ∗ ψ j ( 2 ∇ ψ i ∇ ψ j ∗ ) ] ∑ k ψ k ∗ ψ k = 1 2 ∑ i ∑ j 4 ψ i ∗ ψ j ∇ ψ i ∇ ψ j ∗ / ∑ k ψ k ∗ ψ k , \begin{aligned}
&\frac{\sum_i\sum_j \left [ \psi_i^* \psi_j (2 \nabla \psi_i \nabla \psi_j^*) \right ]}
{\sum_k \psi^*_k\psi_k} \\
& = \frac{1}{2} \sum_i \sum_j 4\psi_i^* \psi_j \nabla \psi_i \nabla \psi_j^* /\sum_k \psi^*_k\psi_k,
\end{aligned} ∑ k ψ k ∗ ψ k ∑ i ∑ j [ ψ i ∗ ψ j ( 2∇ ψ i ∇ ψ j ∗ ) ] = 2 1 i ∑ j ∑ 4 ψ i ∗ ψ j ∇ ψ i ∇ ψ j ∗ / k ∑ ψ k ∗ ψ k , assuming real orbitals (not sure if what would happen otherwise, feel free to double check yourself!), we have:
∑ i ∑ j [ ψ i ∗ ψ j ( 2 ∇ ψ i ∇ ψ j ∗ ) ] ∑ k ψ k ∗ ψ k = 1 2 ∑ i ∑ j ( ∇ ψ i ∇ ψ j ψ i ψ j + ∇ ψ i ψ j ψ i ∇ ψ j + ψ i ∇ ψ j ∇ ψ i ψ j + ψ i ∇ ψ j ∇ ψ i ψ j ) / ∑ k ψ k ψ k = 1 2 ∑ i ( ∇ ψ i ψ i + ψ i ∇ ψ i ) ∑ j ( ∇ ψ j ψ j + ψ j ∇ ψ j ) / ∑ k ψ k ψ k = 1 2 ( ∇ ∑ i ψ i ψ i ) 2 ∑ k ψ k ψ k = 1 2 ( ∇ ρ ) 2 ρ . (18) \begin{aligned}
&\frac{\sum_i\sum_j \left [ \psi_i^* \psi_j (2 \nabla \psi_i \nabla \psi_j^*) \right ]}
{\sum_k \psi^*_k\psi_k} \\
&= \frac{1}{2} \sum_i \sum_j ( \nabla \psi_i \nabla \psi_j \psi_i \psi_j +\nabla \psi_i \psi_j \psi_i \nabla \psi_j + \psi_i \nabla \psi_j \nabla \psi_i \psi_j + \psi_i \nabla \psi_j \nabla \psi_i \psi_j) /\sum_k \psi_k\psi_k \\
&= \frac{1}{2} \sum_i (\nabla \psi_i \psi _i + \psi_i \nabla \psi_i) \sum_j (\nabla \psi_j \psi _j + \psi_j \nabla \psi_j) /\sum_k \psi_k\psi_k\\
&= \frac{1}{2} \frac{(\nabla \sum_i \psi_i \psi_i)^2}{\sum_k \psi_k\psi_k}\\
&= \frac{1}{2} \frac{(\nabla \rho)^2}{\rho}. \tag{18}
\end{aligned} ∑ k ψ k ∗ ψ k ∑ i ∑ j [ ψ i ∗ ψ j ( 2∇ ψ i ∇ ψ j ∗ ) ] = 2 1 i ∑ j ∑ ( ∇ ψ i ∇ ψ j ψ i ψ j + ∇ ψ i ψ j ψ i ∇ ψ j + ψ i ∇ ψ j ∇ ψ i ψ j + ψ i ∇ ψ j ∇ ψ i ψ j ) / k ∑ ψ k ψ k = 2 1 i ∑ ( ∇ ψ i ψ i + ψ i ∇ ψ i ) j ∑ ( ∇ ψ j ψ j + ψ j ∇ ψ j ) / k ∑ ψ k ψ k = 2 1 ∑ k ψ k ψ k ( ∇ ∑ i ψ i ψ i ) 2 = 2 1 ρ ( ∇ ρ ) 2 . ( 18 ) Putting all three terms in Eq. 16 back together, we get:
∇ r ⃗ 2 2 ρ x ( r ⃗ 1 , r ⃗ 2 ) ∣ r ⃗ 2 → r ⃗ 1 = ∇ 2 ρ − 2 τ + 1 2 ( ∇ ρ ) 2 ρ . \nabla_{\vec r_2}^2 \rho_{x}\left(\vec r_1, \vec r_2\right)|_{\vec r_2 \to \vec r_1} = \nabla^2 \rho - 2\tau + \frac{1}{2} \frac{(\nabla \rho)^2}{\rho}. ∇ r 2 2 ρ x ( r 1 , r 2 ) ∣ r 2 → r 1 = ∇ 2 ρ − 2 τ + 2 1 ρ ( ∇ ρ ) 2 . or,
∇ r ⃗ 2 2 ρ x ( r ⃗ 1 , r ⃗ 2 ) ∣ r ⃗ 2 → r ⃗ 1 = ∇ s ⃗ 2 ρ x ( r ⃗ 1 , s ⃗ ) ∣ s ⃗ → 0 = ∇ 2 ρ − 2 τ + 1 2 ( ∇ ρ ) 2 ρ . (19) \nabla_{\vec r_2}^2 \rho_{x}\left(\vec r_1, \vec r_2\right)|_{\vec r_2 \to \vec r_1} = \nabla_{\vec s}^2 \rho_{x}\left(\vec r_1, \vec s\right)|_{\vec s \to 0} = \nabla^2 \rho - 2\tau + \frac{1}{2} \frac{(\nabla \rho)^2}{\rho}. \tag{19} ∇ r 2 2 ρ x ( r 1 , r 2 ) ∣ r 2 → r 1 = ∇ s 2 ρ x ( r 1 , s ) ∣ s → 0 = ∇ 2 ρ − 2 τ + 2 1 ρ ( ∇ ρ ) 2 . ( 19 ) Remember that what we really want is the Taylor expansion of $P^{\sigma \sigma}_\text{cond} (\vec r_1, \vec r_2)$. So, instead of $\left\langle \rho_x (\vec r_1, \vec r_2) \right\rangle$ in Eq. 15, we now write:
⟨ P cond σ σ ( r ⃗ 1 , r ⃗ 2 ) ⟩ = ( 1 + 1 6 s 2 ∇ r ⃗ 2 2 + ⋯ ) P cond σ σ ( r ⃗ 1 , r ⃗ 2 ) ∣ r ⃗ 2 → r ⃗ 1 = P cond σ σ ( r ⃗ 1 , r ⃗ 1 ) + 1 6 s 2 ∇ r ⃗ 2 2 P cond σ σ ( r ⃗ 1 , r ⃗ 2 ) + ⋯ . (20) \begin{aligned}
\left\langle P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{2}\right) \right\rangle&=\left.\left(1+\frac{1}{6} s^{2} \nabla_{\vec r_2}^{2}+\cdots\right) P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{2}\right)\right|_{\vec r_2 \to \vec r_1}\\
&=P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{1}\right) + \frac{1}{6} s^{2} \nabla_{\vec r_2}^{2} P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{2}\right) + \cdots.
\end{aligned} \tag{20} ⟨ P cond σσ ( r 1 , r 2 ) ⟩ = ( 1 + 6 1 s 2 ∇ r 2 2 + ⋯ ) P cond σσ ( r 1 , r 2 ) r 2 → r 1 = P cond σσ ( r 1 , r 1 ) + 6 1 s 2 ∇ r 2 2 P cond σσ ( r 1 , r 2 ) + ⋯ . ( 20 ) The first term is:
P cond σ σ ( r ⃗ 1 , r ⃗ 1 ) , P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{1}\right), P cond σσ ( r 1 , r 1 ) , and we see that it equals to $0$ as:
P cond σ σ ( r ⃗ 1 , r ⃗ 1 ) = ρ σ ( r ⃗ 1 ) − ∣ ρ σ ( r ⃗ 1 , r ⃗ 1 ) ∣ 2 ρ σ ( r ⃗ 1 ) = ρ σ ( r ⃗ 1 ) − ρ σ ( r ⃗ 1 ) = 0. \begin{aligned}
P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{1}\right) &= \rho^{\sigma}(\vec r_{1}) - \frac{\left|\rho^{\sigma}(\vec r_{1},\vec r_{1})\right|^2}{\rho^{\sigma}(\vec r_{1})} \\
&= \rho^{\sigma}(\vec r_{1}) - {\rho^{\sigma}(\vec r_{1})} =0.
\end{aligned} P cond σσ ( r 1 , r 1 ) = ρ σ ( r 1 ) − ρ σ ( r 1 ) ∣ ρ σ ( r 1 , r 1 ) ∣ 2 = ρ σ ( r 1 ) − ρ σ ( r 1 ) = 0. The second term can be easily calculated using Eq. 19:
1 6 s 2 ∇ r ⃗ 2 2 P cond σ σ ( r ⃗ 1 , r ⃗ 2 ) ∣ r 2 → r 1 = 1 6 s 2 [ ∇ 2 ρ − ( ∇ 2 ρ − 2 τ + 1 2 ( ∇ ρ ) 2 ρ ) ] = 1 3 ( τ − 1 4 ( ∇ ρ ) 2 ρ ) s 2 , (21) \begin{aligned}
\frac{1}{6} s^{2} \nabla_{\vec r_2}^{2} P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{2}\right) |_{r_2 \to r_1} &= \frac{1}{6} s^{2} \left[ \nabla^2 \rho - \left( \nabla^2 \rho - 2\tau + \frac{1}{2} \frac{(\nabla \rho)^2}{\rho} \right ) \right]\\
&=\frac{1}{3} \left(\tau - \frac{1}{4} \frac{(\nabla \rho)^2}{\rho} \right ) s^{2}, \tag{21}
\end{aligned} 6 1 s 2 ∇ r 2 2 P cond σσ ( r 1 , r 2 ) ∣ r 2 → r 1 = 6 1 s 2 [ ∇ 2 ρ − ( ∇ 2 ρ − 2 τ + 2 1 ρ ( ∇ ρ ) 2 ) ] = 3 1 ( τ − 4 1 ρ ( ∇ ρ ) 2 ) s 2 , ( 21 ) so that the entire (spherically averaged) Taylor expansion becomes (when $\vec r_2 \to \vec r_1$):
⟨ P cond σ σ ( r ⃗ 1 , r ⃗ 2 ) ⟩ ∣ r ⃗ 2 → r ⃗ 1 = 1 3 [ τ ( r ⃗ 1 ) − 1 4 ( ∇ ρ ( r ⃗ 1 ) ) 2 ρ ( r ⃗ 1 ) ] s 2 + ⋯ . (22) \left\langle P^{\sigma\sigma}_\text{cond}\left(\vec r_{1}, \vec r_{2}\right) \right\rangle |_{\vec r_2 \to \vec r_1} = \frac{1}{3} \left[\tau (\vec r_1) - \frac{1}{4} \frac{(\nabla \rho(\vec r_1))^2}{\rho(\vec r_1)} \right] s^{2} + \cdots. \tag{22} ⟨ P cond σσ ( r 1 , r 2 ) ⟩ ∣ r 2 → r 1 = 3 1 [ τ ( r 1 ) − 4 1 ρ ( r 1 ) ( ∇ ρ ( r 1 ) ) 2 ] s 2 + ⋯ . ( 22 ) The Taylor expansion Eq. 22 of the spherically averaged pair probability succinctly conveys electron localization information. The smaller the probability of finding a second like-spin electron near the reference point, the more highly localized is the reference electron. Hence, electron localization is related to the “smallness” of the expression:
D σ ( r ⃗ 1 ) = τ σ ( r ⃗ 1 ) − 1 4 ( ∇ ρ σ ( r ⃗ 1 ) ) 2 ρ σ ( r ⃗ 1 ) , D_{\sigma}(\vec r_{1})=\tau_{\sigma}(\vec r_{1})-\frac{1}{4} \frac{\left(\nabla \rho_{\sigma}(\vec r_{1})\right)^{2}}{\rho_{\sigma}(\vec r_{1})}, D σ ( r 1 ) = τ σ ( r 1 ) − 4 1 ρ σ ( r 1 ) ( ∇ ρ σ ( r 1 ) ) 2 , or,
D σ ( r ⃗ ) = τ σ ( r ⃗ ) − 1 4 ( ∇ ρ σ ( r ⃗ ) ) 2 ρ σ ( r ⃗ ) , (23) D_{\sigma}(\vec r)=\tau_{\sigma}(\vec r)-\frac{1}{4} \frac{\left(\nabla \rho_{\sigma}(\vec r)\right)^{2}}{\rho_{\sigma}(\vec r)}, \tag{23} D σ ( r ) = τ σ ( r ) − 4 1 ρ σ ( r ) ( ∇ ρ σ ( r ) ) 2 , ( 23 ) where $\sigma$ stands for electrons spin.
One can verify that:
$D_{\sigma}$ is still a probability density, hence it’s non-negative.
$D_{\sigma}$ vanishes in the special case of one-electron systems and hence also vanishes in multielectron systems in regions dominated by a single, localized, $\sigma$-spin orbital.
the relationship between electron localization and Eq. 23 is an “inverse” relationship in the sense that high localizability is implied by small $D_{\sigma}$.
$D_{\sigma}$ is not bounded.
Finally, to assign physical meaning to $D$ and to make it a bounded object, the ELF is defined as:
E L F ( r ⃗ ) = ( 1 + χ σ 2 ( r ⃗ ) ) − 1 , (11) \mathrm{ELF}(\vec r)=\left(1+\chi_{\sigma}^{2}(\vec r)\right)^{-1},
\tag{11} ELF ( r ) = ( 1 + χ σ 2 ( r ) ) − 1 , ( 11 ) where,
χ σ = D σ ( r ⃗ ) / D σ 0 ( r ⃗ ) , (12) \chi_{\sigma}=D_{\sigma}(\vec r) / D_{\sigma}^{0}(\vec r),
\tag{12} χ σ = D σ ( r ) / D σ 0 ( r ) , ( 12 ) and,
D σ 0 ( r ⃗ ) = 3 5 ( 6 π 2 ) 2 / 3 ρ σ ( r ⃗ ) 5 / 3 , (13) D_{\sigma}^{0}(\vec r)=\frac{3}{5}\left(6 \pi^{2}\right)^{2 / 3} \rho_{\sigma}(\vec r)^{5 / 3},
\tag{13} D σ 0 ( r ) = 5 3 ( 6 π 2 ) 2/3 ρ σ ( r ) 5/3 , ( 13 ) where $D_{\sigma}^{0}$ corresponds to a uniform electron gas with charge density equal to the local value of $\rho_{\sigma}(r)$.
One thing to note is that $D_0$ also depends on the local value of $\rho$ so that the dependence is canceled out in ELF and it only reflects the degrees of localization of electrons.
Understanding $D$
We have shown that $D$ which is called the Pauli kinetic density is:
D σ = τ σ − 1 4 ( ∇ ρ σ ) 2 ρ σ , D_{\sigma}=\tau_{\sigma}-\frac{1}{4} \frac{\left(\nabla \rho_{\sigma}\right)^{2}}{\rho_{\sigma}}, D σ = τ σ − 4 1 ρ σ ( ∇ ρ σ ) 2 , where:
τ σ = ∑ i N ∣ ∇ ψ i ∣ 2 , \tau_{\sigma}=\sum_i^N |\nabla \psi_i|^2, τ σ = i ∑ N ∣∇ ψ i ∣ 2 , is twice the positively defined kinetic energy density $\tau_\text{Fermion}$.
Now, if we consider electrons to be bosons and no Pauli exclusion principle exist (feels really wired writing this…), we can write the wavefunctions $\psi_i$ as:
ψ i = 1 N ρ , \psi_i = \frac{1}{\sqrt N} \sqrt \rho, ψ i = N 1 ρ , and so the positively defined kinetic energy density $\tau$ is:
τ Boson = 1 2 ∑ N ∣ ∇ 1 N ρ ∣ 2 = 1 2 ∑ N 1 N ∣ ∇ ρ ∣ 2 = 1 2 ( 1 2 ρ − 1 / 2 ⋅ ∇ ρ ) 2 = ( ∇ ρ ) 2 8 ρ . (B1) \begin{aligned}
\tau_{\text{Boson}} &=\frac{1}{2} \sum^N \left| \nabla \frac{1}{\sqrt N} \sqrt \rho \right|^2\\
&=\frac{1}{2} \sum^N \frac{1}{N} |\nabla \sqrt \rho|^2\\
&=\frac{1}{2} \left( \frac{1}{2} \rho^{-1/2} \cdot \nabla \rho \right) ^2\\
&=\frac{(\nabla \rho)^2}{8 \rho}. \tag{B1}
\end{aligned} τ Boson = 2 1 ∑ N ∇ N 1 ρ 2 = 2 1 ∑ N N 1 ∣∇ ρ ∣ 2 = 2 1 ( 2 1 ρ − 1/2 ⋅ ∇ ρ ) 2 = 8 ρ ( ∇ ρ ) 2 . ( B1 ) Because this is the Bosonic version of the positively defined kinetic energy density (some times also called the von-Weizacker kinetic energy density), it’s the lower bound of the Fermion version:
τ Fermion ≥ τ Boson , \tau_{\text{Fermion}} \geq \tau_{\text{Boson}}, τ Fermion ≥ τ Boson , and the difference between them (to be precisely, twice the difference) is the so-called Pauli kinetic energy $D$:
D = 2 ⋅ τ Fermion − 2 ⋅ τ Boson D = 2\cdot \tau_{\text{Fermion}} - 2\cdot \tau_{\text{Boson}} D = 2 ⋅ τ Fermion − 2 ⋅ τ Boson Difference between $\tau$ and $T$
The kinetic energy density of the system $T$ is defined as:
T = − 1 2 ∑ i ψ i ∗ ∇ 2 ψ i T = -\frac{1}{2}\sum_i \psi_i^* \nabla^2 \psi_i T = − 2 1 i ∑ ψ i ∗ ∇ 2 ψ i Now, let’s consider re-writing $T$ as :
T = − 1 2 ∑ i ψ i ∗ ∇ 2 ψ i = 1 2 ∑ i ∇ ψ ∗ ∇ ψ i − 1 2 ∑ i ∇ ψ ∗ ∇ ψ i − 1 2 ∑ i ψ i ∗ ∇ 2 ψ i = 1 2 ∑ i ∇ ψ ∗ ∇ ψ i − [ 1 2 ∑ i ∇ ψ ∗ ∇ ψ i − 1 4 ∑ i ψ i ∗ ∇ 2 ψ i − 1 4 ∑ i ψ i ∇ 2 ψ i ∗ ] = 1 2 ∑ i ∇ ψ ∗ ∇ ψ i − 1 4 [ ∇ 2 ∑ i ψ i ∗ ψ i ] = 1 2 ∑ i ∇ ψ ∗ ∇ ψ i − 1 4 ∇ 2 ρ = 1 2 ∑ i ∣ ∇ ψ i ∣ 2 − 1 4 ∇ 2 ρ . \begin{aligned}
T &= -\frac{1}{2} \sum_i \psi_i^* \nabla^2 \psi_i \\
&=\frac{1}{2} \sum_i \nabla \psi^* \nabla \psi_i - \frac{1}{2} \sum_i \nabla \psi^* \nabla \psi_i -\frac{1}{2} \sum_i \psi_i^* \nabla^2 \psi_i \\
&=\frac{1}{2} \sum_i \nabla \psi^* \nabla \psi_i - \left[ \frac{1}{2} \sum_i \nabla \psi^* \nabla \psi_i -\frac{1}{4} \sum_i \psi_i^* \nabla^2 \psi_i - \frac{1}{4} \sum_i \psi_i \nabla^2 \psi_i^* \right] \\
&=\frac{1}{2} \sum_i \nabla \psi^* \nabla \psi_i - \frac{1}{4} \left[ \nabla^2 \sum_i \psi_i^* \psi_i \right] \\
&=\frac{1}{2} \sum_i \nabla \psi^* \nabla \psi_i - \frac{1}{4} \nabla^2 \rho \\
&=\frac{1}{2} \sum_i |\nabla \psi_i|^2 - \frac{1}{4} \nabla^2 \rho. \\
\end{aligned} T = − 2 1 i ∑ ψ i ∗ ∇ 2 ψ i = 2 1 i ∑ ∇ ψ ∗ ∇ ψ i − 2 1 i ∑ ∇ ψ ∗ ∇ ψ i − 2 1 i ∑ ψ i ∗ ∇ 2 ψ i = 2 1 i ∑ ∇ ψ ∗ ∇ ψ i − [ 2 1 i ∑ ∇ ψ ∗ ∇ ψ i − 4 1 i ∑ ψ i ∗ ∇ 2 ψ i − 4 1 i ∑ ψ i ∇ 2 ψ i ∗ ] = 2 1 i ∑ ∇ ψ ∗ ∇ ψ i − 4 1 [ ∇ 2 i ∑ ψ i ∗ ψ i ] = 2 1 i ∑ ∇ ψ ∗ ∇ ψ i − 4 1 ∇ 2 ρ = 2 1 i ∑ ∣∇ ψ i ∣ 2 − 4 1 ∇ 2 ρ . The first term is what we called “positively defined kinetic energy” and we have an additional second term.
Actually, if we count the second term in to Eq. 23 we see that it cancels out with the Bosonic part,
D = 2 ⋅ T Fermion − 2 ⋅ T Boson = ∑ i ∣ ∇ ψ i ∣ 2 − 1 2 ∇ 2 ρ − [ 1 4 ( ∇ ρ σ ) 2 ρ σ − − 1 2 ∇ 2 ρ ] , \begin{aligned}
D&=2\cdot T_{\text{Fermion}}-2\cdot T_{\text{Boson}}\\
&= \sum_i |\nabla \psi_i|^2 - \frac{1}{2} \nabla^2 \rho - \left[ \frac{1}{4} \frac{\left(\nabla \rho_{\sigma}\right)^{2}}{\rho_{\sigma}} - - \frac{1}{2} \nabla^2 \rho \right],
\end{aligned} D = 2 ⋅ T Fermion − 2 ⋅ T Boson = i ∑ ∣∇ ψ i ∣ 2 − 2 1 ∇ 2 ρ − [ 4 1 ρ σ ( ∇ ρ σ ) 2 − − 2 1 ∇ 2 ρ ] , where we have used Eq. B1 to get $T_{\text{Boson}}$.
Free electron gas $D_0$
Assuming periodic boundry condition, in 3D, the Schordinger Equation has solution:
ψ k ⃗ ( r ⃗ ) = e i k ⃗ r ⃗ , (A1) \psi_{\vec k}(\vec r) = e^{i \vec k \vec r}, \tag{A1} ψ k ( r ) = e i k r , ( A1 ) where $\vec k$ takes the following values:
k x = 2 π L n x , k y = 2 π L n y , k z = 2 π L n z , n x , n y , n z = 0 , ± 1 , ± 2 , ⋯ . (A2) \begin{aligned}
k_x =\frac{2\pi}{L} n_x, \quad k_y &=\frac{2\pi}{L} n_y, \quad k_z =\frac{2\pi}{L} n_z, \\
n_x, n_y, n_z &= 0, \pm 1,\pm 2, \cdots. \tag{A2}
\end{aligned} k x = L 2 π n x , k y n x , n y , n z = L 2 π n y , k z = L 2 π n z , = 0 , ± 1 , ± 2 , ⋯ . ( A2 ) The dispersion relation is:
ε = ℏ 2 2 m ⋅ k ⃗ 2 , (A3) \varepsilon=\frac{\hbar^2}{2m} \cdot \vec{k}^2, \tag{A3} ε = 2 m ℏ 2 ⋅ k 2 , ( A3 ) or in atomic unit:
ε = 1 2 k ⃗ 2 . (A4) \varepsilon=\frac{1}{2} \vec{k}^2. \tag{A4} ε = 2 1 k 2 . ( A4 ) Since we are considering Fermi gas, each state (only differ by their $k$ value) can only take one electron (for one spin channel). Starting from the lowest energy state ($\vec k=0$) the highest filled state has the Fermi vector $\vec k_F$.
From Eq. A2, we see that each kpoint takes a volume of $(\frac{2\pi}{L})^3$. The volume of occupied “k-sphere” is $\int_0^{k_F} dk = \frac{4}{3} \pi k_F^3$, If we only consider one spin channel, the total number of allowed occupied states is:
N = 4 3 π k F 3 ( 2 π L ) 3 = 4 3 π k F 3 ( 2 π ) 3 V , (A5) N = \frac{\frac{4}{3}\pi k_F^3}{(\frac{2\pi}{L})^3}= \frac{\frac{4}{3}\pi k_F^3}{\frac{(2\pi)^3}{V}}, \tag{A5} N = ( L 2 π ) 3 3 4 π k F 3 = V ( 2 π ) 3 3 4 π k F 3 , ( A5 ) where $V=L^3$ is the volume of the unitcell. Note that the number of occupied states equals to the number of electrons in that spin channel $N=N_e$.
From Eq. A4 we can calculate the Fermi vector $k_F$:
k F = ( 6 π 2 N e V ) 1 3 . (A6) k_F = \left( 6 \pi^2 \frac{N_e}{V} \right)^{\frac{1}{3}}. \tag{A6} k F = ( 6 π 2 V N e ) 3 1 . ( A6 ) To calculate the average kinetic energy per electron(or total energy per electron since there are no potential energy term in the free electron SE) , we need to sum up energy of each occupied state and then devided that by the number of electrons in the system,
T = 1 N e ∑ k ε k = ( ∫ 0 k ⃗ F d k ⃗ ( 2 π ) 3 / V ) − 1 1 ( 2 π ) 3 / V ∫ 0 k ⃗ F 1 2 k 2 d k ⃗ = ( 4 π 3 k F 3 ) − 1 ⋅ 2 π 5 k F 5 = 1 2 ⋅ 3 5 ( 6 π 2 N e V ) 2 3 . (A7) \begin{aligned}
T&=\frac{1}{N_e} \sum_k \varepsilon_k \\
&=\left(\frac{\int_0^{\vec k_F} d\vec k}{(2\pi)^3/V}\right)^{-1} \frac{1}{(2\pi)^3/V} \int_0^{\vec k_F} \frac{1}{2}k^2 d\vec k\\
&=\left( \frac{4\pi}{3} k_F^3 \right)^{-1} \cdot \frac{2\pi}{5} k_F^5\\
&=\frac{1}{2} \cdot \frac{3}{5} \left( 6 \pi^2 \frac{N_e}{V} \right)^{\frac{2}{3}}. \tag{A7}
\end{aligned} T = N e 1 k ∑ ε k = ( 2 π ) 3 / V ∫ 0 k F d k − 1 ( 2 π ) 3 / V 1 ∫ 0 k F 2 1 k 2 d k = ( 3 4 π k F 3 ) − 1 ⋅ 5 2 π k F 5 = 2 1 ⋅ 5 3 ( 6 π 2 V N e ) 3 2 . ( A7 ) If we define the local charge density as $\rho=N_e/V$, Eq. A7 becomes:
T = 1 2 ⋅ 3 5 ( 6 π 2 ) 2 3 ( ρ ) 2 3 . T = \frac{1}{2} \cdot \frac{3}{5} \left( 6 \pi^2 \right)^{\frac{2}{3}} \left( \rho \right)^{\frac{2}{3}}. T = 2 1 ⋅ 5 3 ( 6 π 2 ) 3 2 ( ρ ) 3 2 . Since $T$ is the kinetic energy per electron, the kinetic energy at specific point can be calculated by:
T ρ ( r ⃗ ) = 1 2 ⋅ 3 5 ( 6 π 2 ) 2 3 ( ρ ) 5 3 . (A8) T \rho(\vec r) = \frac{1}{2} \cdot \frac{3}{5} \left( 6 \pi^2 \right)^{\frac{2}{3}} \left( \rho \right)^{\frac{5}{3}}. \tag{A8} Tρ ( r ) = 2 1 ⋅ 5 3 ( 6 π 2 ) 3 2 ( ρ ) 3 5 . ( A8 ) Finally, since:
τ σ = ∣ ∇ ψ i ∣ 2 \tau_{\sigma} =|\nabla \psi_i|^2 τ σ = ∣∇ ψ i ∣ 2 is actually twice of the kinetic energy density, we need to double the value we just calculated to get $\tau_{\sigma}$:
τ σ = 3 5 ( 6 π 2 ) 2 3 ( ρ ) 5 3 (A9) \tau_{\sigma} = \frac{3}{5} \left( 6 \pi^2 \right)^{\frac{2}{3}} \left( \rho \right)^{\frac{5}{3}} \tag{A9} τ σ = 5 3 ( 6 π 2 ) 3 2 ( ρ ) 3 5 ( A9 )
REF0,Eq.1;