View source for Electronic Properties of Materials/Quantum Mechanics for Engineers/Many Electron Atoms

This is the tenth chapter of the first section of the book Electronic Properties of Materials.



We now have a solutions for Hydrogen.<math display="block">\begin{align}
H &= T_{electron} + V_{ion-electric}
\\ &= {- \hbar^2 \over 2m} \nabla^2 - {|e|^2 \over r}
\\ H \psi_{nlmm_s} (r) &= E_n \psi_{nlmm_s}
\\ E_n &= {-1 \over n^2} \reals
\\ \psi_{nlmm_s} &= R_{nl} Y_{lm} \chi_{m_s} 
\end{align}</math> 


In undergraduate chemistry class, we were shown this solution and told "... and so it goes for the rest of the periodic table..." However, the truth is, it isn't so simple. 

Consider Lithium:

<FIGURE> "Diagram of Lithium Ion" (Description)

We can apply center of mass corrections but for simplicity, assume that the nucleus has <math>M_{noc} = \infty</math>.

What does hydrogen look like?

<math display="block">H = T_{electron} + V_{ion-electron} + V_{electron-electron}</math><math display="block">\begin{align}
T_{electron} &= T_{e_1} + T_{e_2} + T_{e_3}
\\ &= {-\hbar^2 \over 2m} \left ( 
    {\nabla_1}^2 + {\nabla_2}^2 + {\nabla_3}^2 \right )
\\
\\ V_{ion-electron} &= {-3e^2 \over r_1} + {-3e^2 \over r_2} + {-3e^2 \over r_3}
\\ V_{electron-electron} &= {e^2 \over |r_1 - r_2|} + {e^2 \over |r_2 - r_3|} 
    + {e^2 \over |r_3 - r_1|}
\end{align}</math>


Without the electron-electron interaction the problem would be much simpler.<math display="block">\begin{align}
H &= T_1 + V_{ie1} + T_2 + V_{ie2} + T_3 + V_{ie3}
\\ &= h_1 + h_2 + h_3
\end{align}</math>

This makes for a separable solution: 

<math>\begin{align}
H\Psi &= E\Psi
\\ \Psi &= \psi_1 \psi_2 \psi_3
\end{align}</math>

Here, each <math>\psi_i
</math> terms is just the hydrogen wave function with a slight modification for <math>V= {3 e^2 \over r} = {z e^2 \over r}</math>.

<math>\begin{align}
    E = \varepsilon_1 + \varepsilon_2 + \varepsilon_3
    \\ \varepsilon_i = -{Z^2 \mathbb{R} \over \eta^2}
\end{align}</math>

Here <math>Z</math> is the atomic number. 

Since electrons are fermions, <math>\Psi</math> must be totally anti-symmetric as seen in 2 Particles in a Box in Chapter 4.

<math display="block">\Psi(r_1, r_2, r_3) = {1 \over \sqrt{6}} \left | 
    \begin{align}
        \psi_1(r_1) \quad \psi_2(r_1) \quad \psi_3(r_1)
        \\ \psi_1(r_2) \quad \psi_2(r_2) \quad \psi_3(r_2)
        \\ \psi_1(r_3) \quad \psi_2(r_3) \quad \psi_3(r_3)
    \end{align}
\right |</math>

This is known as the Slater Determinant. 

=== Approximating Lithium ===
Unfortunately, we can't just discard electron to electron interactions. The magnitude of <math display="inline">V_{ion-electron}</math> and <math display="inline">V_{electron-electron}</math> are equivalent. So, what do we do? Solve the problem as is? Well, after you solve <math>H_{Li}</math>, you could quite possibly earn a Nobel Prize and/or a Fields Medal which makes this a wholly impractical method for this course. Instead we will use approximate methods for quantum mechanical problems more complicated than our earlier Hydrogen solution, and we currently have two primary methods for doing this.

===== Thomas-Fermi Model =====
The Thomas-Fermi Model is a semi-classical statistical method where we replace the exact potential <math display="inline">V_{ion-electron}</math> and <math display="inline">V_{electron-electron}</math> with an effective potential, <math display="inline">V_{eff}</math>, the screened coulomb potential. 

Consider an electron gas with so many electrons that you don't have to count them one at a time, but instead just consider a continuous charge distribution. If we put a positive charge in the electron cloud then the positive charge would attract the electrons so the electron distribution would no longer be uniform. That said, the negative charges nearest the positive charge "screen" the electrons further out such that the fringe electrons don't feel as strong of an attraction. 

Working out the details of this method requires a discussion of electron gasses, provided later in this text. Meanwhile, you should know that this hypothetical situation turns into a spherically distributed electron gas where the charge density, <math>\rho</math>, is slowly varying. Although <math>\rho(r)</math> does change with <math>r</math>, at each <math>r</math>, we can treat <math>\rho</math> as a uniform gas. Additionally, a point charge <math>+Q</math>, placed in a uniform electron gas has a screened potential of:

<math display="block">\Phi(r)= {Q \over r} e^{-k_or}</math>

This is sometimes called a Yukawa Potential. Also, <math display="inline">{1 \over k_o}</math> is the Thomas-Fermi screening length and is related to the density of quantum states of the highest filled state, called the density of states at the Fermi Level. 

Finally, for charge neutral atoms, where the number of electrons equals the atomic number, with a slowly varying electron gas in a spherically symetric central potential:

<math display="block">V(r) \simeq \underbrace{-{Ze^2 \over r}}_\text{Attractive Potential} + 
    \underbrace{1.794 \ {Z^{4 \over 3} \over a_o}}_\text{Repulsive Potential} + \ \cdots</math>

<FIGURE> "Electron Gas Potential" (Description)

The truth of the matter is that the Thomas-Fermi Model isn't a great solution for the multi-electron atom. Back in the 1920's it was quite successful and is still useful for simple approximations and the initial input to more advanced methods. 

===== Hartree-Fock Method =====
A better method is the iterative Hartree-Fock Method, or the Self-Consistent Field Method. For this method, we start by assuming that we can write:

<math display="block">\Psi = \prod^n_{i=1} \psi_i</math>This assumption further implies that the can say that <math display="inline">H = \sum^N_{i=1} h_i</math>. So, what is a reasonable <math display="inline">h_i</math>?

<math display="block">h_i = {-\hbar^2 \over 2m} \nabla^2_i + V_o (r_i)+{1 \over 2} \sum_{j \neq i} {e^2 \over |r_i - r_j|}</math>

Now, this gives us:

<math display="block">H = \sum^N_{i=1} \left [ 
       {-\hbar^2 \over 2m} \nabla^2_i + V_o (r_i)+{1 \over 2} \sum_{j \neq i} {e^2 \over |r_i - r_j|}
\right ]</math>

<math display="inline">V_o(r_i)</math> is all the potential terms '''not''' involving electron-electron interaction. For the single atom, <math display="inline">V_o(r_i) = {Ze^2 \over r_i}</math>, but this method is generally applicable to more complex systems such as molecules. Hartree-Fock is the intellectual forbearer of several modern quantum chemical techniques. 

<nowiki><Source> "Primer on Calculus of Variation and Lagrange Multipliers"</nowiki>

=== Energy Function ===
The functional we are interested in is the energy, where <math display="inline">E \left [ \Psi(r) \right ] = \int \Psi^* H \Psi \ dr</math>. We want to find the extrema subject to the constraints <math display="inline">\int \psi^*_i \psi_i \ dr =1</math> for all <math>i</math>. We will do this by using Lagrangian multipliers <math>\varepsilon_i</math> such that the integral we're interested in is:

<math display="block">\int \Psi^* H \Psi \ \operatorname{d}\!r 
    - \sum_i \varepsilon_i \int \psi^*_i \psi_i \ \operatorname{d}\!r</math>

==== Applying Calculus of Variation ====
We want to use calculus of variation to find the functions that make this function stationary (an extrema). Consider this function:<math display="block">I[\phi_1(x) \ \phi_2(x)] = \int f(\phi_1(x), \ \phi_2(x), \ x) \ \operatorname{d}\!x </math>We want to vary the function to look for the stability condition. In other words, we want to change <math>\phi_i(x) </math> at an arbitrary constant <math>x </math>.<math display="block">\Delta I = I[\phi_1(x) + \varepsilon_1 \eta_2 (x), \ \phi_2(x) + \varepsilon_2 \eta_2 (x)]
    - I[\phi_1, \ \phi_2] </math>

Here <math>\varepsilon_i </math> is a small constant (we take the limit of <math>\varepsilon_i </math> as it goes to zero), and <math>\eta_i </math> is an arbitrary function that deviates <math>\phi_i </math> (this must be continuous over the limits of integration. This gives us:<math display="block">\Delta I = 
    \int f(\phi_1(x) + \varepsilon_1 \eta_2 (x), \ \phi_2(x) + \varepsilon_2 \eta_2 (x), \ x] - 
    \int f(\phi_1(x), \ \phi_2(x), \ x) </math>

The Maclaurin expansion in powers of <math>\varepsilon_i </math> is:<math display="block">f(x) = f(0) + f'(0) + {1 \over 2!}f''(0) + {1 \over 3!} f'''(0) + \cdots </math><math display="block">\begin{align}
\Delta I &= 
    \int f(\phi_1(x), \ \phi_2(x), \ x) + 
        {\partial f \over \partial \phi_1} \ \varepsilon_1 \eta_1 (x) + 
        {\partial f \over \partial \phi_2} \ \varepsilon_2 \eta_2 (x) + 
        \text{Higher Order Terms} - f(\phi_1(x), \ \phi_2(x), x) 
\\ &= {\partial f \over \partial \phi_1} \ \varepsilon_1 \eta_1 (x) + 
        {\partial f \over \partial \phi_2} \ \varepsilon_2 \eta_2 (x)
\\ \delta I &= 
    \int {\partial f \over \partial \phi_1} \ \delta \phi_1 +
        {\partial f \over \partial \phi_2} \ \delta \phi_2
\\ &= \int \delta f \qquad \text{(The variation of the functional} \ f )
\end{align} </math>

<math>\delta \phi_i = \varepsilon_i \eta_i  </math> in the limit as <math>\varepsilon_i </math> goes to zero.


Let's return to our problem:

<math>I = \int \Psi^* H \Psi - \sum \varepsilon_j \psi^*_j \psi_j \operatorname{d}\!r </math>

Writing this function out in detail gives us:<math display="block">I = \int \underbrace{(\psi^*_1, \ \cdots, \ \psi^*_n)(h_1, \ \cdots, \ h_n)(\psi_1, \ \cdots, \ \psi_n) -
    \varepsilon_1 \psi^*_1 \psi_1 - \cdots - \varepsilon_n \psi^*_n \psi_n}_{f(\psi^*_1, \ \cdots, \psi^*_n)} \operatorname{d}\!r </math>

Note that you could alternatively carry through <math>\psi </math> in this step, but in the end you don't need them. I chose to utilize <math>\psi^* </math> instead of <math>\psi </math> for purely aesthetic reasons. This gives us:

<math display="block">\delta I = 0 = \int \sum^n_{q=1} {\partial f \over \partial \psi^*_q} \delta \psi^*_q \ \operatorname{d}\!r </math>

Looking at the <math>q^{th} </math> term in the summation, <math display="inline">{\partial f \over \partial \psi_q^*} \ \delta \psi_q^* </math>, we can expand to:

<math display="block">\begin{align}
    & \sum_{j=1}^n \left [ 
        (\psi_1^*, \ \cdots, \ \delta \psi^*_q, \ \cdots, \psi_n^*)(h_j)
        (\psi_1, \ \cdots, \ \psi_q, \ \cdots, \psi_n) \right ] - \delta \psi_q^* \varepsilon_q \psi_q 
    \\ & \sum_{j=1}^n \left [ 
        (\psi_1^*, \ \cdots, \ \delta \psi^*_q, \ \cdots, \psi_n^*) 
            \left[ \underbrace{{-\hbar^2 \over 2m} \nabla_j^2}_{T_j} + \underbrace{V_o (r_j)}_{V_j^o} 
                + \underbrace{{1 \over 2} \sum_{k \neq j}^n {e^2 \over |r_j - r_k|}}_{V_{jk}^{ee}} \right ] 
            (\psi_1, \ \cdots, \ \psi_q, \ \cdots, \psi_n)\right ] - \delta \psi_q^* \varepsilon_q \psi_q 
\end{align} </math>

For the <math>T_j </math> and <math>V_j^o </math> terms, there exists two conditions:

# When <math>j \neq q </math>, the terms inside the summation become <math>math </math>, but integrating this will go to zero because <math>\delta \psi_q^* </math> and <math>\psi_q </math> are orthogonal.
# The only non-zero terms come about for <math>j = q </math>

Thus, dividing through by <math>\Psi^* \Psi </math> gives us <math>\delta \psi_q^* (T_q + V_q^o) \psi_q </math>. Here, the <math>V_{jk}^{ee} </math> term is again zero for <math>j \neq q </math>, because the orthogonal functions come out from <math>\delta \psi_q^* \psi_q </math>. The resulting term is:

<math display="block">\delta \psi_q^* \left ( 
    {1 \over 2} \sum _{k \neq q} \psi_k^* {e^2 \over |r_q - r_k|} \psi_k \right ) \psi_q </math>

Because <math>V_{jk}^{ee} </math> depends on both <math>j_q </math> and <math>k </math>, we retain the wavefunctions <math>\psi_k </math> inside the summation.

The functional variation now looks like:

<math display="block">\begin{align}
\delta I &= \int \sum_{q=1}^n \delta \psi_q^* (T_q + V_q^o) \psi_q + \delta \psi_q^* 
    \left ( {1 \over 2} \sum_{k \neq q} \psi_k^* {e^2 \over |r_q - r_k|} \psi_k \right ) \psi_q
    - \delta \psi_q^* \varepsilon_q \psi_q
\\ &=0
\end{align} </math>

Which can be solved as <math>n </math> equations.

<math display="block">\left [ {- \hbar^2 \over 2m} \nabla^2 + V_o (r) + V_d(r) \right ] \psi_k(r) 
    = \varepsilon_k \psi_k(r) </math>

Where <math>V_d </math> is the Hartree Term which is:    <math>V_d = \sum_{j \neq k}^n \int \psi_j^* (r_j) \psi_j(r_j) {e^2 \over |r_j - r_k|} \ \operatorname{d}\!r </math>

This looks like a Schrodinger equation, but is it? Yes, sort of... <math>\psi_i </math> are Hartree wave functions. Remember that these are single electron <math>\psi_i </math> since we started with <math>H = \sum_i h_i </math>. <math>\varepsilon_i  </math> are Lagrangian multipliers for enforce normalization. However, we can use these <math>\{ \psi_i \} </math> and <math>\varepsilon_i  </math> to study the system.

===== Hartree Equation =====
The Hartree equation has a physically intuitive form (if anything in quantum mechanics is physically intuitive). THe first term is kinetic energy, the second term is the electron-ion or any external potential energy (such as E or B), and the third term is the electron-electron potential energy defined by the sum:

<math display="block">\sum_{j \neq k} \int \underbrace{\psi^*_j(r_j) \psi_j(r_j)}_{\rho_j} 
    {e^2 \over |r_j - r_k|} \ \operatorname{d}\!r</math>

Here <math>\rho_j</math> is the charge density of the <math>j^{th}</math> electron. So this sum can be thought of as an electron at <math>r_k</math> interacting with some charge at <math>r_j</math>. 

What's wrong with this picture? Well, excluding our single particle assumption, we're putting several electrons into this system, which are fermions, but our wavefunction is not totally anti-symmetric. What we need to do is write <math>\Psi</math> using a Slater determinant. 

<math display="block">\begin{align}
\Psi &= {1 \over \sqrt{n!}} \left | \begin{matrix} 
    \psi_1(r_1) & \cdots & \psi_n(r_1)
    \\ \cdots
    \\ \psi_1(r_n) & \cdots & \psi_n(r_n)
\end{matrix} \right |
\\ &= {1 \over \sqrt{n!}} \sum_{\rho} (-1)^\rho \psi_1(r_1) \ \cdots \ \psi_n(r_n)
\end{align}</math>

Put this back into our original energy function and apply calculus of variation methods to the whole expression to find our answer. This was done by Fock (and Slater) in 1930. The resulting set of expressions is exactly the same except for one additional term. <math display="block">\left [ {-\hbar^2 \over 2m} \nabla^2 + V_o(r) + V_d(r) + V_{\text{exchange}}(r) \right ]
    \psi_k(r) = \varepsilon_k \psi_k(r)</math>

The evaluation of <math>V_{\text{exchange}}</math> is quite tricky. The expression for the expectation value is:<math display="block">\langle E_{\text{exchange}} \rangle = {1 \over 2}
    \sum_{1\neq j} \int \int \psi_i^*(r_i) \psi_j^*(r_j) {e^2 \over |r_i - r_j|}
    \underbrace{\psi_i(r_j) \psi_j(r_i)}_{\operatorname{d}\!{r_i} \operatorname{d}\!{r_j}} </math>

In many cases it is acceptable to approximate using the so-called "free electron exchange" <math display="inline">V_{\text{exchange}}(r) = -e^2 \left ( {3 \over \pi} \right )^{1 \over 3} \rho(r)^{1 \over 3}</math>

===== Self-Consistent Field Loop =====
These equations are not trivial to solve. Notice that the operator depends on <math>\rho(r)</math>, but to get <math>\rho(r)</math> you need <math>\psi_i(r)</math> which requires solving the set of equations. This is where we solve by iterative, self-consistent methods. This technique is the basis for all of modern quantum mechanical methods. 

Write <math display="inline">(T + V(r)) \psi_k(r_k) = \varepsilon_k \psi_k(r_k)</math>, where <math display="inline">V(r)</math> is all the terms, not <math>T</math>. It will depend on <math>\rho(r)</math> or <math>\psi(r)</math> or both. 

# Guess <math display="inline">V^{(1)}(r)</math>. You can use known solutions, random guesses, the Thomas-Fermi, or any other method of guessing you wish.
# Solve the system of equations to find <math>\{\psi_i^{(1)}(r)\}</math>.
# Calculate  from <math>\{\psi_i^{(1)}(r)\}</math>.
# If <math display="inline">V^{(2)} \simeq V^{(1)}</math> end. Otherwise return to step 2 and keep looping until <math display="inline">V^{(n)} \simeq V^{(n+1)}</math> to within a given tolerance.

This is called a self-consistent field loop (S.C.F.).

These type of SCF methods can be applied to all sorts of systems, but since this chapter is about many-electron problems, we will look at the Herman-Skillman atomic data, calculated from the Hartree-Fock method in 1963.

<nowiki><SOURCE> "Herman-Skillman" atomic data"</nowiki>

{{navigation|Book=Electronic Properties of Materials
|current=Many Electron Atoms
|previous=Quantum Mechanics for Engineers/Perturbation Methods
|next=Quantum Mechanics for Engineers/Density Functional Theory}}{{BookCat}}