5. Beyond hydrogen: the electronic structure of atoms, molecules, and solids#
While we can solve for the eigenstates and corresponding energies of the electron in the hydrogen atom exactly, we’re out of our depth for essentially any other system[1]. To see why other systems are so difficult, let’s consider the two simplest cases: the helium atom, and the hydrogen molecule ion. The helium atom consists of a single nucleus with charge \(+2e\), and two electrons with charge \(-e\). The interaction between the nucleus and each of the electrons is much the same as in hydrogen, acting radially as well. The interaction between the two electrons however depends on their separation, and thus adds a non-spherically symmetric term to the Hamiltonian. Likewise, the hydrogen molecule ion (an ionized hydrogen molecule) consists of three particles: two single-proton nuclei with charge \(+e\) and one electron with charge \(-e\). As the protons are much heavier than the electron, we can work in the approximation that they are stationary (known as the Born-Oppenheimer approximation), but even then, we have to account for the fact that both nuclei interact with the electron. As they can’t both be at the origin, we again get a Hamiltonian that is not spherically symmetric. While this specific case can still be solved by a coordinate transformation in which the Hamiltonian separates, the math becomes increasingly difficult, and no solutions are known for more complicated systems.
Although attempts at solving for the wavefunctions (or even only the energies) of non-spherically symmetric Hamiltonians are doomed[2], we can get quite close by employing various approximations. In addition to the Born-Oppenheimer one introduced above, the main approach involves the variational principle, which allows you to construct an upper bound to the ground state energy of any Hamiltonian, and in some cases to the energies of exited states as well.
5.1. The variational principle#
Theorem 5.1
For any Hamiltonian \(\hat{H}\), the ground state energy \(E_1\) is bounded from above by
where \(\psi(x)\) is an arbitrary function in the Hilbert space \(L^2(\infty)\).
Proof. Finding a state with energy larger than \(E_1\) is trivial (by construction, any of the excited states have higher energy). To show that no state has lower energy, recall that the eigenstates \(\psi_n(x)\) and corresponding eigenvalues \(E_n\) of the Hamiltonian (the solutions of \(\hat{H} \psi_n = E_n \psi_n\)) are complete and orthonormal, so \(\braket{\psi_m | \psi_n} = \delta_{mn}\), and for any \(\psi \in L^2(\infty)\) we can write
Moreover, from the normalization of \(\psi(x)\) we have
We now use these properties to estimate the expectation value of \(\hat{H}\) in the state \(\psi(x)\):
which is exactly equation (5.1).
The variational principle has an extremely useful consequence: to estimate the energy of the ground state, you can use any guess for the wavefunction \(\psi(x)\), and it will always result in an upper bound. Therefore, by optimizing your guess, i.e., tuning it so that the expectation value of the Hamiltonian comes out as low as possible, you can construct approximations of the actual ground state energy.
Example 5.1 (an estimate for the ground state energy of the harmonic oscillator)
To see how the variational principle works, let’s apply it to a simple example: the one-dimensional harmonic oscillator Hamiltonian:
The first step is to select a good trial function. What we want is a function for which we can actually calculate the expectation value of the energy, so the integrals of this function should be fairly easy. Moreover, we want the function to have a parameter that we can vary, so we can find the lowest possible estimate of the expectation value of the Hamiltonian. A popular choice for such a trial function is a Gaussian, by its virtue of having very easy integrals. We take
where \(b\) is our parameter, and \(A = (2b/\pi)^{1/4}\) is the normalization constant. Splitting the Hamiltonian into a kinetic and a potential part, we get
so
Now as the variational principle states that the ground state energy of \(\hat{H}\) is less than or equal to this estimate for any value of the parameter \(b\), we can minimize the expression in equation (5.2) to get the lowest possible upper bound. We have
Substituting this value for \(b\) into our estimate, we find:
Equation (5.3) happens to give the exact result, which is because the actual ground state of the harmonic oscillator happens to be a Gaussian (with the found optimal value of \(b\)). That of course is just a lucky guess, though it illustrates the fact that the variational principle can get you close to the actual value, without ever having to solve a differential equation.
5.1.1. First application: the ground state energy of the helium atom#
The Hamiltonian of the helium atom contains two kinetic terms (one for each electron), two identical potential terms for the interaction between each electron and the nucleus, and a single potential term for the interaction between the two electrons. The terms for the kinetic plus nuclear interaction energy of each electron are the same as those for hydrogen (except with a \(+2e\) nuclear charge), so we can write:
where the numbers 1 and 2 refer to the two electrons. As stated above, it’s the interaction term between the two electrons that’s the source of all trouble; if it were not there, the solution would simply be the product of the hydrogen wave functions for each of the electrons. In particular, the ground state would be
where we’ve corrected for the \(+2e\) charge of the helium nucleus. The corresponding energy of each of the individual electrons would be \(4E_1\), where \(E_1\) is the ground state energy of hydrogen; for two electrons this would add up to \(8E_1 = -109\;\mathrm{eV}\) (see Exercise 5.1a). The actual (experimentally determined) value is \(-79\;\mathrm{eV}\), so the estimate is far off. In fact, it underestimates the actual value. We can easily do better, by applying the variational principle: we can use the ‘product ground state’ of equation (5.5) as our guess and calculate the expectation value of the actual Hamiltonian, given by equation (5.4)b. Again, the only difficulty is the electron-electron interaction part (see Exercise 5.1); we now find an estimate above (and much closer to) the experimental value.
As the variational principle allows us to pick any wave function, we could probably do better than just the product of the hydrogen ground state wave functions. A well-known phenomenon in electrostatics is screening: if in between yourself and a large positive charge there is a negative charge, the effective charge you feel is less than the total positive charge. Suppose that somehow the two electrons screen the charge of the nucleus to each other; they’d then feel an effective nuclear charge \(+Ze\) (where, presumably, \(Z<2\)). Let’s try that and optimize for \(Z\). We re-write the (total, unchanged) Hamiltonian as
where the first two parts now combine to ‘hydrogen atoms’ with nuclear charge \(Ze\), and hence ground-state energies \(Z^2 E_1\). The corresponding hydrogen atom ground state product eigenstate is closely related to equation (5.5):
The expectation value of the helium Hamiltonian in this state can be calculated from values we’ve computed before:
Optimizing for \(Z\) gives \(Z=\frac{27}{16} \approx 1.69\), which gives us an estimate of \(E_1 \leq \Braket{\hat{H}} = -77.5\;\mathrm{eV}\). Much better than just the product state, and, once we’ve done the hard integrals (which we fortunately have to do only once), not much more work.
5.2. Molecules#
5.2.1. Covalent bonds#
Covalent bonds are the things that hold molecules together. No doubt you’ve encountered them in chemistry as the ‘links’ between atoms, with different atoms having different numbers of available link slots. For example, with its four link slots (and being the lightest, and therefore most abundant, such type of atom), carbon forms the backbone of all of organic chemistry, including the amino acids that together form the proteins that make life work. A carbon atom can link to four other atoms, as it does in methane (\(\mathrm{CH}_4\)), but also form double or even triple bonds with other multi-valent atoms, for example making two double bonds with two oxygen atoms in carbon dioxide (\(\mathrm{CO}_2\)). As we’ll see in this section, covalent bonds are a purely quantum mechanical phenomenon. We will first work out the general concept in one dimension; in the sections that follow, we’ll see some concrete examples in 3D.
For a one-dimensional system, suppose we have two quantum particles that can both be in single-particle states \(\psi_a(x)\) and \(\psi_b(x)\); we’ll take these states to be normalized and orthogonal. If these were classical particles, we could paint one red and one blue, and thus distinguish them, allowing us to create a combined state \(\psi(x_1, x_2) = \psi_a(x_1) \psi_b(x_2)\). For indistinguishable quantum-mechanical particles, this wavefunction will not do, as probabilities should be invariant if we swap the particles; we therefore end up with either a symmetric or an antisymmetric combination (see Section 3.6):
Suppose now that our two particles represent valence electrons in two atoms that are relatively close together. The nuclei of these atoms are heavy compared to the electrons; we’ll take their positions fixed for now (we’ll use the separation as a parameter in a variational calculation later); we’ll also consider other non-valence electrons to be ‘strongly bound’ to their nucleus. We can then calculate the mean distance between the two electrons, or (more practically), the mean of the square of their separation. If the electrons would be distinguishable particles, the calculation of this expectation value is straightforward, and gives
where the subscripts \(a\) and \(b\) indicate the state that the particle is in (i.e., \(\Braket{x^2}_a = \braket{\psi_a(x)|x^2|\psi_a(x)}\), etc.). For indistinguishable particles, the calculation is not very hard either, but gets us more terms:
where we have a new ‘exchange term’, in which the two single-particle wavefunctions both contribute
Due to the exchange term, identical particles will, on average, be found either closer together or further apart than distinguishable ones. We get an attractive ‘exchange force’ for bosons (with a plus sign in equation (5.9), but also for electrons in an antisymmetric spin singlet (and thus symmetric orbital) configuration. These forces originate from the extra exchange term which is purely quantum-mechanical. By itself, the force does of course not guarantee that the particles will be bound together; that will only happen if the resulting energy is also lower than for the configuration where they are apart. Therefore, rather than calculating separations, we should be calculating energies. We’ll do so in the sections below, and we’ll find similar quantum-mechanical terms to the exchange term above; it is these terms that cause the reduction of energy of ‘joint states’ compared to separated states, and thus are responsible for making covalent bonds stable.
5.2.2. The hydrogen-molecule ion and the LCAO#
The simplest two-atomic molecule is hydrogen (\(\mathrm{H}_2\)): two single-proton nuclei with two electrons in shared orbits. Molecular hydrogen can be ionized, creating \(\mathrm{H}_2^{+}\): a molecular ion with only one electron. As the protons are much heavier than the electron (by a factor of about \(2000\)), to good approximation, we can take them as stationary, and consider only the motion of the electron. This approach is known as the Born-Oppenheimer approximation. For a given distance between the nuclei, we can solve for the position of the single electron exactly, because its Hamiltonian is separable when expressed in the right (here elliptical, see Section 5.4) coordinates. This exact solution is however significantly more complicated than the (admittedly already complicated) solutions for the hydrogen atom, and as we’re typically more interested in the neutral form (with two electrons), not particularly useful. Alternatively, we can approximate the solutions of the hydrogen molecule ion using a technique that is widely used in molecular chemistry, and which will, using the variational principle again, allow us to calculate an estimate of the separation of the two protons. This new approximation is known as the LCAO: the linear combination of atomic orbitals.
Fig. 5.1 Coordinate system for the molecular hydrogen ion. The two nuclei are located at \(z=\pm \frac12 R\). The single electron (red dot) is a distance \(\bm{r}_1\) from the left nucleus, and a distance \(\bm{r}_2\) from the right nucleus. The easiest coordinate system to work in for this configuration is an elliptical one, as described in Section 5.4.#
For our molecular hydrogen ion, the Hamiltonian is given by
where we’ve taken the protons fixed at \(z = \pm \frac12 R\) and \(r_i\) denotes the distance of the electron to the \(i\)th proton, see Fig. 5.1. Note that if we did not fix the protons, we’d have to add kinetic terms for each of them, so the Born-Oppenheimer approximation has simplified our problem significantly. Moreover, we’ll treat the separation \(R\) between the protons as the parameter in our variational approach, so for calculating the wave function of the electron, the last term of (5.11) simply adds a constant to the total energy, equal to
where, as before, \(E_1\) is the ground state energy of the hydrogen atom, and \(a\) the Bohr radius.
Instead of solving for the eigenfunctions of the Hamiltonian in equation (5.11) exactly, we’ll use the variational method to estimate the ground state energy of our system. Our trial function is now going to be a linear combination of (hydrogen) atomic ground states[3], one for each of the nuclei:
Note that \(\bm{r}_1\) now indicates the distance to the first nucleus, and \(\bm{r}_2\) that to the second; both \(\bm{r}_1\) and \(\bm{r}_2\) are functions of the general position \(\bm{r}\). The key difference with the trial function (5.5) for the helium atom is that instead of a product, we now have a sum of atomic orbitals. In some sense, that makes life easier, but as we’ll see shortly, we still get a bunch of non-trivial cross terms. The first cross terms already appear if we simply try to normalize our trial function; we have
The first two terms in equation (5.13) give us no trouble (the atomic ground state functions are already normalized), but the third term is nontrivial. It occurs frequently enough that it got its own name: the overlap integral, as it measures how much the wave functions \(\psi_{100}(\bm{r}_1)\) and \(\psi_{100}(\bm{r}_2)\) overlap in space. To evaluate it, we go to elliptical coordinates (see Section 5.4), in which the integral becomes easy. We have[4]
Given the value of the overlap integral \(S\), the normalization constant \(A\) in equation (5.12) becomes \(A = 1/\sqrt{2(1+S)}\).
The second quantity to calculate is the expectation value of the Hamiltionian, equation (5.11), in the trial state (5.12). As the Hamiltonian consists of the atomic hydrogen Hamiltonian plus an additional term for the second nucleus and a constant for the potential between the nuclei, we have
and similarly for \(\psi_{100}(\bm{r}_2)\). The calculation of the expectation value of \(\hat{H}\) will therefore give us some more cross terms:
where we used that \(E_1 = - (e^2/4\pi\varepsilon_0) \cdot (1/2a)\), with \(a\) the Bohr radius. Note that we absorbed the Bohr radius in the definition of \(D\) and \(X\) to make them dimensionless.
In equation (5.15) we introduced two new integrals[5] the direct integral \(D\) and the exchange integral \(X\). The direct integral emerges from the interaction between the (charged) electron in an atomic orbit centered on the first nucleus and the presence of the (charged) second nucleus. The exchange integral comes about due to the interaction of the overlap of the charges of the two (atomic) orbitals and one of the nuclei. Unlike the direct and overlap integral, which would also appear classically, the exchange integral has no classical analog, as classically the electron cannot be in a state that is a superposition of two orbitals of the two nuclei.
Calculating the direct and exchange integrals is straightforward using elliptical coordinates, and proceeds in a very similar fashion to the calculation of the overlap integral in equation (5.14), see Exercise 5.2. The result is given by
In terms of the various integrals, we now find an easy expression for \(\braket{\hat{H}}\):
The difference between the total energy and that of the hydrogen atom ground state is plotted in Fig. 5.2. The figure also shows the energy we get if we use the other obvious linear combination of atomic orbitals, which is antisymmetric in the two nuclei (indicated with a minus subscript; we’ll use the plus subscript for the symmetric combination in equation (5.12)):
The calculation of the expectation value of the Hamiltonian in the antisymmetric state goes completely analogous to the symmetric state (see problem Section 5.\ref{prob:hydrogenionintegrals}); the result is
As it turns out, we have \(D > X > 0\) for the molecular hydrogen ion. Consequently, the antisymmetric state always has an energy larger than \(E_1\), while the symmetric state has a region where its energy is below \(E_1\), with a minimum at a separation of approximately \(2.5 a\). The hydrogen molecule ion is thus stable[6] it will not spontaneously fall apart into a hydrogen atom and a free proton, as the symmetric state including the two protons has a lower total energy.
Show code cell source
%config InlineBackend.figure_formats = ['svg']
import numpy as np
import matplotlib.pyplot as plt
from myst_nb import glue
def S(s):
return (1 + s + s*s/3) * np.exp(-s)
def D(s):
return (1/s) - (1+(1/s)) * np.exp(-2*s)
def X(s):
return (1+s) * np.exp(-s)
def Eplus(s):
return (2/s) - 2 * (D(s) + X(s)) / (1 + S(s))
def Eminus(s):
return (2/s) - 2 * (D(s) - X(s)) / (1 - S(s))
# Plot energies
s = np.linspace(.1, 5, 500)
fig, ax = plt.subplots(figsize=(6,4))
line1 = ax.plot(s, Eplus(s), label='$E_{+}$')
line2 = ax.plot(s, Eminus(s), label='$E_{-}$')
ax.axhline(y = 0, color = 'k', linestyle = ':')
ax.set_xlim(0,5)
ax.set_ylim(-0.2,1.5)
ax.set_xlabel('$R/a_0$', fontsize=16)
ax.set_ylabel('$\\frac{E-E_1}{E_1}$', fontsize=16)
ax.legend(fontsize=16)
# Save graph to load in figure later (special Jupyter Book feature)
glue("hydrogenmoleculeionLCAOenergy", fig, display=False)
Fig. 5.2 Energy of the LCAO approximation for the hydrogen molecule ion for the symmetric \(\psi_+\) (blue line) and antisymmetric \(\psi_-\) (orange line) combination of the two atomic orbitals, as a function of the separation of the nuclei (\(2R\)), in units of the Bohr radius \(a\). For the symmetric case, we find a state with an energy lower than the energy of an electron in a hydrogen atom, so the hydrogen molecule ion is stable.#
5.2.3. Variational principle for general molecules#
Inspired by the success of the variational principle in predicting the existence of a stable molecular hydrogen ion, we may use the same LCAO approximation to study other molecules. We’ll assume that only the valence electrons (the electrons in the highest occupied but not completely filled ‘shell’, where a ‘shell’ consists of orbitals with the same principal quantum number \(n\)) participate. As our trial function, we’ll then write general linear combinations of all possible atomic orbitals available to our valence electron(s). We thus write our trial function as
where the \(\psi_i\) are atomic orbitals for one of the constituent molecules. As in Section 5.2.2, we’ll encounter two types of integrals when evaluating the expectation value of the Hamiltonian in our trial function. The first type is due to the overlap between various atomic orbitals, generalizing the overlap integral \(S\) of equation (5.14). The second type stems from direct and cross terms in the expectation value of the Hamiltonian itself, generalizing the direct and exchange integrals of equations (5.16) and (5.17). As we’ll encounter such terms for every pair of constituent atomic wavefunctions \(\psi_i\), \(\psi_k\), the corresponding integrals are conveniently arranged in a matrix and known as the overlap matrix elements, given by
and the Hamiltonian matrix elements, defined as
We can easily express the energy of the trial function (5.21) in terms of these matrix elements:
By the variational principle, the energy of the actual ground state of the Hamiltonian[7] is less than this estimate. Consequently, we can freely tune the coefficients \(c_i\) in our trial function (5.21) to minimize the energy in equation (5.24). At the minimum, the variation of the energy must vanish, which gives[8]
which implies that \(E\) is minimized under the condition that
Equation (5.26) is a matrix equation: it tells us that the solutions (collections \(\lbrace c_k \rbrace\) of coefficients) span the nullspace of the matrix \(H_{ik} - E S_{ik}\), or in other words, that the matrix \(H_{ik} - E S_{ik}\) is not invertible[9]. For such solutions to exist, the determinant of that matrix must vanish, so a necessary condition for the energy to be minimized is that
5.2.4. The hydrogen molecule#
We can use the general principle derived in Section 5.2.3 to re-derive the result for the hydrogen molecule ion (see Exercise 5.2). More interestingly, we can also use it to optimize a trial function for the energy of the hydrogen molecule. To get the Hamiltonian for that molecule, we add a kinetic term for a second electron, and a repulsive term between the two electrons, to the Hamiltonian of the hydrogen molecule ion:
In equation (5.28), the numbers label the electrons, and the letters the nuclei; thus \(r_1^A\) is the distance between electron 1 and nucleus A, see Fig. 5.3. Like for the helium atom, for the most part, the Hamiltonian of the hydrogen molecule is simply twice that of the hydrogen molecule ion; the only additional (but nasty) term comes from the interaction between the two electrons.
Fig. 5.3 Coordinate system for molecular hydrogen, with two protons labeled with letters and two electrons labeled with numbers.#
We now have multiple ways of constructing a trial function for the hydrogen molecule. The perhaps conceptually simplest one is to look for states that are product states of the two underlying hydrogen atoms. We can make linear combinations of these as we did for the molecular hydrogen ion; these combinations could in principle be symmetric (with a singlet spin state for overall antisymmetry) or antisymmetric (with a triplet spin state). These combinations are known as the Heitler-London approximation, with a sign to indicate whether the combination is symmetric or antisymmetric:
where \(\sigma_{2}(1, 2)\) represents the singlet spin state and \(\sigma_{+}(1, 2)\) the triplet spin state. Finding the normalization constant for this combination is straightforward, and results in
where \(S\) is the overlap integral we defined in equation (5.14). The expectation value of the energy again contains the direct and exchange integrals, but also two new integrals, which we get from calculating the contribution of the electron-election interaction to the energy. We label these integrals \(J\) and \(K\); they are defined in equations (5.33)a and (5.33)b below. The calculation of the expectation value of the Hamiltonian in the Heitler-London approximation is somewhat tedious but otherwise straightforward; the result is given by:
where \(E_1\) is still the ground state energy of hydrogen, and we used the relation \(e^2/24 \pi \varepsilon_0 = - 2 a E_1\). The two energies are plotted in Fig. 5.4. As we can read off from the plot, the symmetric orbital combination with singlet spin state has an energy lower than the two separate atoms, and thus predicts that molecular hydrogen will be stable. It however overestimates the distance between the nuclei (\(1.64\) times the Bohr radius, where the experimentally obtained value is \(1.40\)) and underestimates the binding energy (the minimum is at about \(-0.23 E_1\), while the experimental value is \(-0.35 E_1\); experimental values from [Huber and Herzberg, 1979], page 250).
Instead of starting from two atoms, we could also start from our best approximation to the ground state of the molecular hydrogen ion. As we found in Section 5.2.2 that for a single electron the symmetric orbital has the lowest energy, a natural starting point would be to put both electrons in such a symmetric orbital. Fortunately, we can do so without violating the Pauli exclusion principle, as long as we put the two electrons in the singlet spin state. Our trial function (which I’ll call the two-ion approximation) then reads
where \(\psi_{+}(\bm{r}_i)\) means that we take the symmetric combination of atomic orbitals for electron \(i\) (as we did in equation (5.12) for the molecular hydrogen ion) and \(\sigma_{-}(1, 2)\) again represents the singlet spin state. As we already know the expectation value of the hydrogen molecule ion in the \(\psi_{+}\) state, and the only additional term in the Hamiltonian is the interaction between the two electrons, the energy of our trial state will be twice that of the hydrogen molecule ion plus a correction term. Because \(\psi_{+}(\bm{r})\) is itself a linear combination of two terms, and we now have two electrons, the correction term will have four subterms, all contributing to the calculation of \(\Braket{1/|\bm{r}_1 - \bm{r}_2|}\):
where \(r_{12} = |\bm{r}_1 - \bm{r}_2|\) and we again included the Bohr radius \(a\) in the definition to make our integrals dimensionless[10]. Equation (5.33)a represents the electron-electron repulsion and (5.33)b, (5.33)c and (5.33)d the various overlap charge repulsions. Evaluating the integrals in equation (5.33) is feasible, though painful (the full expressions are given in Section 5.5). In terms of these integrals (and those from Section 5.2.2), the energy of the molecule is given by
The resulting curve as a function of nuclear separation is plotted in Fig. 5.4. As you’ll no doubt immediately spot, something must be wrong, as the expectation value of the energy is larger than that of two separate atoms if we push the nuclei to infinity. The reason for this behavior is that our chosen trial function cannot be written as the product of two orbital wave functions, as we started from the solutions of the molecular hydrogen ion. That does not make this solution useless though: by Theorem 5.1, as long as our trial function is in the Hilbert space, the expectation value we get for the Hamiltonian is guaranteed to be larger than the actual energy of the ground state. Unfortunately, the smallest value we get from this two-ion approximation is \(-0.20 E_1\), so the guess is less good than that of the Heitler-London approximation.
In equation (5.32) we chose the symmetric combination of atomic orbital wavefunctions for each of the two electrons (forcing us to put them in the singlet state to make the total state antisymmetric). Like for the molecular hydrogen ion, we could have made other choices: we could have put one of the electrons in the antisymmetric orbital, or both electrons in the antisymmetric orbitals, creating an overall symmetric orbital which again requires an antisymmetric (singlet) spin. If we have two or more states with the same symmetry (we’ll call these configurations, combinations of orbital and spin states), the system could exist in a linear combination of these states, which could be a better approximation of ground state of the whole molecule. Chemists refer to this refinement of the estimate of the ground state energy as configuration interaction: the presence of multiple configurations with the same symmetry affects the energy of the lowest-energy configuration, usually lowering it somewhat.
To distinguish between the various states, chemists have come up with an elaborate naming scheme. We already encountered the notation for atomic orbitals of the form \((nl)^k\), where \(n\) (a number) is the principal quantum number, \(l\) (denoted as a letter, with \(s=0\), \(p=1\), \(d=2\) and \(f=3\)) is the orbital quantum number, and \(k\) the number of electrons in that state. For orbitals in a molecule, the notation is similar, replacing the Roman with Greek letters (\(s \to \sigma\), \(p \to \pi\), etc.), dropping the parentheses, and adding a subscript that indicates whether the orbital configuration[11] is symmetric (or even, indicated with a g, from the German ‘gerade’) or antisymmetric (or odd, indicated with a u, from the German ‘ungerade’). The configuration in equation (5.32) is then written as \(1 \sigma_\mathrm{g}^2\). If we replace one of the symmetric LCAO orbitals with the asymmetric one, we get the configuration \(1 \sigma_\mathrm{g} 1 \sigma_\mathrm{u}\), and making them both asymmetric gives \(1 \sigma_\mathrm{u}^2\). For the whole molecule, they also generalize the \({}^{2S+1}L_J\) notation, keeping the \(S\) and \(L\) as they are (while again replacing Roman by Greek letters) but replacing the \(J\) for the total angular momentum with the (orbital) symmetry label (i.e., \(\mathrm{g}\) or \(\mathrm{u}\)). Therefore, \({}^1\Sigma_\mathrm{g}\) means that we have a wavefunction with total spin \(0\) (the superscript \(1\) following from \(2 \cdot 0 + 1\)), total orbital angular momentum \(0\) (as indicated by the capital \(\Sigma\)), and an even combination of atomic orbitals (meaning that if we swap them, they don’t flip sign), as indicated by the \(\mathrm{g}\) subscript. To further illustrate, the four options we get for the hydrogen molecule are
For the second and third option, we put the electrons in a symmetric (option 2) and antisymmetric (option 3) combination of orbitals, so they should be in either the singlet (5.35)b or the triplet (5.35)c spin state; in neither case the total configuration has the same symmetry as the (supposed) ground state (5.35)a. The ‘doubly antisymmetric’ orbital configuration (5.35)d does have the same symmetry as the ground state, and thus could be put in a linear combination with it. Its energy can be expressed in the same integrals as that of the doubly symmetric combination, and reads
The ground state wavefunction including the configuration interaction then reads
where as before the capital letters label the nuclei and the number indexes the electrons. Because we can tune the coefficients \(c_1\) and \(c_4\) in equation (5.37), we can lower our estimate of the ground state energy, see Fig. 5.4. The lowest energy is given by the combination
The estimate of the ground-state energy we get from this configuration interaction wavefunction does (marginally) better than the Heitler-London approximation: \(-0.24 E_1\).
We have one more card left to play: like for the helium atom, we may expect that the two electrons screen the charge of the two protons. Indeed they do, and repeating our calculations with the screening effect as a free parameter allows us to further refine the estimate of both the ground state energy (to \(-0.30 E_1\)) and the separation between the nuclei (to \(1.43\) times the Bohr radius) [Halpern and Glendening, 2013]. While this technique thus gives a significant improvement, all approximations we’ve tried result in the same conclusion: molecular hydrogen is stable.
Show code cell source
%config InlineBackend.figure_formats = ['svg']
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sc
from myst_nb import glue
# Helper functions / integrals
def S(s):
return (1 + s + s*s/3) * np.exp(-s)
def D(s):
return (1/s) - (1+(1/s)) * np.exp(-2*s)
def X(s):
return (1+s) * np.exp(-s)
def J(s):
return (1/s) - (1/s + 11/8 + 3*s/4 + s*s/6) * np.exp(-2*s)
def A(s):
return (6/s) * ((np.euler_gamma + np.log(s)) * S(s) * S(s) + sc.expi(-4*s) * S(-s) * S(-s) - 2 * sc.expi(-2 * s) * S(s) * S(-s))
def B(s):
return (-25/8 + 23*s/4 + 3 * s * s + s*s*s/3) * np.exp(-2*s)
def K(s):
return (A(s) - B(s))/5
def M(s):
return 5/8
def L(s):
return (s + 0.125 + 5 / (16*s)) * np.exp(-s) - (0.125 + 5 / (16*s)) * np.exp(-3*s)
# Energy estimates
def EHLplus(s):
# Energy from Heitler-London approximation.
# We return the energy minus 2 E1 and with an overall minus sign as the value of E1 is negative.
return (2 / s) - 2 * (2 * D(s) - J(s) + 2 * X(s) * S(s) - K(s)) / (1 + S(s) * S(s))
def EHLminus(s):
return (2 / s) - 2 * (2 * D(s) - J(s) - 2 * X(s) * S(s) + K(s)) / (1 - S(s) * S(s))
def Etwoionplus(s):
# Energy from taking the product of two hydrogen molecule ion states.
return (2/s) - 4 * (D(s) + X(s)) / (1+S(s)) + (J(s) + 2 * K(s) + M(s) + 4 * L(s)) / ((1+S(s)) * (1+S(s)))
def Etwoionminus(s):
return (2/s) - 4 * (D(s) - X(s)) / (1-S(s)) + (J(s) + 2 * K(s) + M(s) - 4 * L(s)) / ((1-S(s)) * (1-S(s)))
def H12(s):
# Additional term in the configuration interaction expectation value.
return (1/((1-S(s)*S(s)))) * (M(s) - J(s))
def ECIplus(s):
# Energy from configuration interaction optimization.
return 0.5* (Etwoionplus(s) + Etwoionminus(s) - np.sqrt((Etwoionplus(s) - Etwoionminus(s))*(Etwoionplus(s) - Etwoionminus(s)) + 4 * H12(s) * H12(s)))
def ECIminus(s):
return 0.5* (Etwoionplus(s) + Etwoionminus(s) + np.sqrt((Etwoionplus(s) - Etwoionminus(s))*(Etwoionplus(s) - Etwoionminus(s)) + 4 * H12(s) * H12(s)))
# Comparison of five energy estimates.
fig,ax = plt.subplots(figsize=(10,8))
lw = 2
fs2 = 14
s = np.linspace(0.1, 5, 500)
ax.plot(s, EHLplus(s), lw = lw, label=r'$E_{+}^{\mathrm{HL}}$')
ax.plot(s, EHLminus(s), lw = lw, label=r'$E_{-}^{\mathrm{HL}}$')
ax.plot(s, Etwoionplus(s), lw = lw, label=r'$E_{+}^{\mathrm{2ion}}$')
ax.plot(s, Etwoionminus(s), lw = lw, label=r'$E_{-}^{\mathrm{2ion}}$')
ax.plot(s, ECIplus(s), lw = lw, label=r'$E^{\mathrm{CI}}$')
ax.set_ylim(-0.25, 1)
ax.set_xlim(0,5)
ax.set_xlabel(r'$R/a$', fontsize = fs2)
ax.set_ylabel(r'$\frac{E - 2 E_1}{E_1}$', fontsize = fs2)
ax.legend(fontsize = fs2)
hline = {"color":"k","linestyle":":","linewidth":lw}
ax.axhline(y = 0, **hline)
# Save graph to load in figure later (special Jupyter Book feature)
glue("hydrogenmoleculeLCAO", fig, display=False)
Fig. 5.4 Various estimates of the energy of the hydrogen molecule: the Heitler-London approximation (equation (5.31)), the combination of two molecular hydrogen ion states (equations (5.34) and (5.36)) and the configuration interaction estimate (equation (5.38)). The energies are plotted as a function of the separation \(R\) between the two nuclei, in units of the Bohr radius \(a\). The energies are the difference between the estimated energy and the energy of two separate hydrogen atoms, in units of the hydrogen atom ground state energy.#
5.2.5. Diatomic molecules#
For homonuclear diatomic molecules (i.e., molecules consisting of two identical atoms) the procedure for constructing estimates for the ground state energy goes analogous to that of the hydrogen molecule. Similar to what we’ve seen for the periodic table, the valence electrons are typically the relevant ones; electrons in completely filled shells are typically closely bound to a nucleus and have little overlap with the electrons of other atoms in a molecule. Moreover, if the energies between two orbitals (centered on the two constituent atoms) are very different, their interaction will have little effect on the overall energy levels of the molecule. To see why that is the case, we consider the general minimization condition on the energy from the variational principle, equation (5.27). Taking \(H_{11} = \alpha_\mathrm{A}\) and \(H_{22} = \alpha_\mathrm{B}\) for the self-interactions, \(S_{12} = S_{21} = S\) for the overlap and \(H_{12} = H_{21} = \beta\), the determinant becomes
Now if \(S = 0\) and \(|\alpha_\mathrm{A} - \alpha_\mathrm{B}| \gg \beta\), the solutions of equation (5.39) are
so they equal the atomic orbital energies (\(\alpha_\mathrm{A}\) and \(\alpha_\mathrm{B}\)) with small corrections \(\pm \beta^2/(\alpha_\mathrm{B}-\alpha_\mathrm{A})\). The strongest bonds (largest energy difference with atomic orbitals) will therefore be between orbitals with comparable energy and significant overlap. Electrons in such orbitals are known as valence electrons. To calculate the energy of a bound state, we should therefore take interactions between all orbitals in a shell (i.e., same \(n\) quantum number) with the same symmetry into account. For example, for the atoms in the second row of the periodic table, the corresponding diatomic molecules will have orbitals that are rotationally symmetric about the axis through the two atoms. These orbitals (known as \(\sigma\) orbitals to chemists) will be linear combinations of the \(\psi_{200}\) and \(\psi_{210}\) atomic orbitals for each of the nuclei. Likewise, there can be orbitals that carry a net angular momentum about the central axis (called \(\pi\) orbitals for \(L=1\)), which are linear combinations of the \(\psi_{21\pm1}\) atomic orbitals.
For each of the eight elements in the second row of the periodic table (the \(n=2\) elements, lithium to neon), we can construct the molecular orbitals by combining the 2s (\(n=2, l=0\)) and the 2p (\(n=2, l=1\)) orbitals. Unsurprisingly, the lowest-energy molecular orbitals are the combinations of two 2s atomic orbitals, in even and odd combinations. At this point, unfortunately, the naming conventions start to deviate. Some authors simply start numbering the orbitals by the order in which they occur in the molecule, making the symmetric and antisymmetric combinations of the 2s orbitals the \(1 \sigma_\mathrm{g}\) and \(1 \sigma_\mathrm{u}\) orbitals. The argument is then that the electrons in the 1s (\(n=1, l=0\)) orbitals don’t contribute to the molecular orbitals. Others however reserve the \(1 \sigma\) orbitals for electrons with \(n=1\), and call the combinations of the 2s orbitals \(2 \sigma_\mathrm{g}\) and \(2 \sigma_\mathrm{u}\). As this choice seems to be the most prevalent, we will follow it here. For \(l=1\), we can make two more linear combinations with \(m_l = 0\), called \(3 \sigma_\mathrm{g}\) and \(3 \sigma_\mathrm{u}\), and four linear combinations of \(m_l = \pm 1\) atomic orbitals. These turn out to be pairwise degenerate, with the antisymmetric combination having lower energy than the symmetric one (because of the inherent antisymmetry of the atomic orbitals around the axis through the two nuclei). These orbitals are known as \(1\pi_\mathrm{u}\) and \(1 \pi_\mathrm{g}\). Fig. 5.5 shows all eight options, with the various diatomic molecules indicated at their highest-occupied molecular orbital (often abbreviated as HOMO). Note that the order of the energy levels differs for the first five and last three elements.
From Fig. 5.5, we can read off the electron configurations in each of the diatomic molecules. For example, in dinitrogen with 10 valence electrons, we occupy the states \((2\sigma_\mathrm{g})^2 (2\sigma_\mathrm{u})^2 (1\pi_\mathrm{u})^4 (3\sigma_\mathrm{g})^2\). All orbital states are occupied by two electrons in a spin singlet, and the total state of the molecule is \({}^1\Sigma_\mathrm{g}\). In contrast, for dioxygen with 12 valence electrons, we read off that we have \((2\sigma_\mathrm{g})^2 (2\sigma_\mathrm{u})^2 (3\sigma_\mathrm{g})^2 (1\pi_\mathrm{u})^4 (1\pi_\mathrm{g})^2\). The top level, \(1\pi_\mathrm{g}\), is now not fully occupied, as we only use two out of the four available orbitals. The electrons can therefore exist in separate orbitals, and thus can have parallel (triplet) spin. It turns out that that configuration has a lower energy than the identical orbital-singlet spin one, so the overall state of the oxygen molecule is \({}^3\Sigma_\mathrm{g}\).
Fig. 5.5 The molecular orbital energy levels of diatomic molecules of the elements in the second row of the periodic table (\(n=2\)). The labels indicate the highest occupied level of the given molecule. The left column indicates the energy levels up to di-nitrogen, the right column of di-oxygen and above.#
5.2.6. Polyatomic molecules#
For polyatomic molecules, naively combining the various available valence electrons into an LCAO gets out of hand quickly. For example, for water, we have two hydrogen atoms, each providing a single electron, and an oxygen, providing four valence electrons. The full LCAO thus has six electrons, leading to a \(6 \times 6\) secular determinant in equation (5.27). Fortunately, with a clever choice of basis, this determinant will factorize, making life significantly easier.
As a standard basis for the LCAO, let us take \(\psi_{100}(r_i)\) (\(i=1, 2\)) for each of the hydrogens (with \(r_i\) the distance to the \(i\)th nucleus) and the set \(\psi_{200}(r_3), \psi_{210}(r_3), \psi_{211}(r_3), \psi_{21-1}(r_3)\) for the oxygen (with distance \(r_3\) to the nucleus). The trick is now to find linear combinations of wavefunctions with the same symmetry (known as symmetry adapted linear combination, SALC), as wavefunctions with different symmetry are perpendicular, and thus their cross terms will drop out of the determinant. For the water molecule, we have an obvious mirror symmetry, so a natural starting point is the hydrogen combination
where A and B are the two hydrogen atoms. There are two atomic orbitals of the oxygen atom that have the same mirror symmetry: \(\psi_{200}(\bm{r}_3)\) and \(\psi_{210}(\bm{r}_3)\). We can also make an antisymmetric combination of the hydrogen atoms:
which has the same symmetry as the \(2p_y\) combination of the oxygen orbitals, given by
Finally, the \(2p_x\) combination of the oxygen orbitals, given by
has no partner of the same symmetry in the hydrogens. Taking as our basis now the combination
the determinant factorizes into a \(3 \times 3\), a \(2 \times 2\) and a \(1 \times 1\) block, and is easily evaluated.
Fig. 5.6 Water molecule. (a) Symmetries of the orbitals of the valence electrons in the oxygen atom (top four pictures) and the two hydrogen atoms (bottom two pictures). Regions with the same blue color are symmetric, regions where blue and orange are juxtaposed are antisymmetric. (b) Resulting energy levels. In the ground state, the lowest four levels are occupied, each by two electrons.#
5.2.7. The Hückel approximation for large molecules#
Even with symmetry-adapted linear combinations, large molecules still tend to generate large determinants, especially when we have large numbers of identical atoms. Therefore, we sometimes need to make further simplifying assumptions. A common one that works well for the \(\pi\)-orbitals of a carbohydrate chain or ring is the Hückel approximation, where we set
All overlap integrals to zero, so \(S_{ik} = \delta_{ik}\) (this is a rather poor approximation for nearest neighbors and the largest source of error in this method).
All diagonal elements of the Hamiltonian to the same value, \(H_{kk} = \alpha\) (quite decent for long carbon chains).
All off-diagonal elements to zero, except for nearest-neighbor terms, which are all set to the same value \(\beta\).
For example, for the four carbons in (1,3)-butadiene (\(\mathrm{H}_2\mathrm{C} = \mathrm{CH} - \mathrm{CH} = \mathrm{CH}_2\)), we get
which gives (for \(x = \alpha-E\)):
For the energies we then obtain
The energies and corresponding wavefunctions are plotted in Fig. 5.7. Perhaps unsurprisingly, we find that the four energy levels correspond to four wavelengths, essentially the first four modes of a wave that ‘fits’ on the ‘string’ made by the four carbon atoms.
Fig. 5.7 Energies and corresponding wavefunctions of the p-orbitals of butadiene.#
5.3. Solids#
Molecules are formed through covalent bonds between atoms. Pictorially, we can say that the atoms share their valence electrons, thus filling up their outermost shell. This description works well for the elements on the right-hand side of the periodic table. The elements in the rightmost column already have full shells, and thus do not form any bonds; they instead form mono-atomic ‘noble’ gases. The elements in the columns just left of the noble gases form one to four covalent bonds with other elements from the same ‘environment’, giving rise to a very rich variety of chemical compounds. Especially relevant are the carbon atoms, which can form four bonds, and constitute the backbone of most polymers, and thus, of life.
For the metals on the left-hand side of the periodic table, the picture is somewhat different. Rather than forming covalent bonds, at (relatively) low temperatures they preferentially arrange in a crystal lattice, with the nuclei and most of the associated electrons highly localized on grid points. Such lattices can consist of a single type of atom (pure metals), multiple different metal atoms (alloys), or combinations of ions of a metal and a non-metal (salts). For now, we’ll focus on the simplest case, of a crystal containing a single type of atom at well-spaced intervals. The valence electrons of these atoms are not just shared with the nearest neighbors; in an extension of what we’ve seen for longer carbon chains in Section 5.2.7, the electrons roam (more or less) freely through the entire solid. Consequently, their states will be highly de-localized, and we should describe them from the perspective of the solid as a whole, rather than as the sum of individual atoms.
5.3.1. The free electron gas model#
The simplest model we can construct for a solid is to ignore the interactions between the valence electrons and the nuclei altogether. If we also ignore the electrostatic interactions between the various electrons, they become, in effect, particles of an ideal gas. As no doubt you know from classical thermodynamics, a classical ideal gas can be described with a handful of parameters (pressure, temperature, number of particles). The energy of the gas in a box can be calculated by considering the only interaction that remains, the collisions the gas molecules make with the boundaries of the box. In the free electron gas model of Sommerfeld we do the same thing for our (quantum-mechanical) electrons: we only consider the interactions with the boundaries of the solid. As the electrons cannot escape the solid, and have no interactions inside it, we thus return to the very simplest case: we treat the solid as an infinite well. The difference with Section 2.2.1 is that we will now work with a three-dimensional well, and that our electrons still have to obey the Pauli exclusion principle.
For a three-dimensional infinite well with size \(L_x \times L_y \times L_z\), solving the time-independent Schrödinger equation is easy, as the equation separates in three one-dimensional infinite well equations for the \(x\), \(y\) and \(z\) components. The solutions to the Schrödinger equation are then simply the products of the three-one dimensional solutions, while the energies of the solutions are the sums of the respective one-dimensional energies. In particular, we write
\end{subequations}
where we have defined the wave vector \(\bm{k} = (k_x, k_y, k_z)\) and \(k = |\bm{k}|\), which, as before, is directly related to the electron’s momentum through \(\bm{p} = \hbar \bm{k}\). Each wave vector corresponds to a point in k-space (also known as reciprocal space or momentum space), the dual space of the coordinate space spanned by the position basis vectors \(\hat{x}\), \(\hat{y}\) and \(\hat{z}\). In this dual space, the points are \(\pi / L\) apart (where \(L\) can be \(L_x\), \(L_y\), or \(L_z\); each point thus has an associated volume
where \(V\) is the (real-space) volume of the solid. Unsurprisingly, we find that the larger our solid is, the closer together the k-space points are; after all, the wave vector is proportional to the inverse of the wave length, and in a larger solid, we can fit larger waves.
Fig. 5.8 Two dimensional representation of the points in k-space corresponding to the various states of the electrons in the free electron gas model. Each points stands for a unique wavevector \(\bm{k}\). The distance between two points is \(\pi/L\), so the volume associated to each point is \(\pi^3/V\). As the energy is proportional to the length of the wavevector squared, points with the same distance to the origin have the same energy. Due to the Pauli exclusion principle, for \(N\) atoms with \(q\) valence electrons each, \(Nq/2\) states will be occupied in the ground state. The surface bounding the volume representing all occupied states in the ground state is known as the Fermi surface; in the free electron gas model, this surface is a piece of a sphere with radius \(k_\mathrm{F}\). To find the total energy in the ground state, we integrate over shells with thickness \(\mathrm{d}k\), in which all states have the same energy.#
The electrons of our free electron gas are the valence electrons of the constituent atoms. Let \(q\) be the number of valence electrons per atom, and \(N\) be the number of atoms. While \(q\) is usually small, \(N\) is (very) large, on the order of Avogadro’s number for a macroscopic solid. If electrons were bosons, the ground state of the free electron gas would be the one in which all of them are in \(\psi_{111}\). However, because electrons are spin-\(\frac12\) fermions, we can only put two electrons in that state. Above the ground state, there are three degenerate states with the same energy, \(\psi_{211}\), \(\psi_{121}\), and \(\psi_{112}\), also three with the next energy state, and so on. To determine the order of the states, we can look at k-space: the energy is proportional to the square of the wave vector, and thus states with equal energy have the same distance to the origin in k-space. To find the ground-state energy of our free electron gas, we put two electrons in each block in k-space, going radially outward, until we’ve placed all our electron. The ground-state energy is then the sum of all the energies of all the blocks. The shape of the collection of filled sites in k-space will be an octant of a solid sphere; for large numbers of particles, the approximation of the shape as a sphere will be highly accurate (see Fig. 5.8 for a two-dimensional representation). Finding the highest occupied energy level is then easy: we simply equate the volume of the octant of the sphere to the total volume of all occupied sites:
where \(\rho = Nq/V\) is the free electron density (the number of valence electrons per unit volume), and \(k_\mathrm{F}\) is the radius of the Fermi sphere. All models for solids have a ground state in which the electrons occupy a volume in k-space; the boundary of that volume is known as the Fermi surface; our sphere is just the simplest example. Calculating the energy of the highest-energy electrons in the ground state (the Fermi energy) is now trivial:
The Fermi energy is the energy of the most energetic electrons. To find the total energy of all electrons in the ground state we need to sum over all occupied sites in k-space. With the large numbers we’re dealing with, we can approximate this sum by an integral. As we’ve seen, all states with equal distance to the origin in k-space have the same energy; therefore, we only need to integrate over shells (thin spheres) with radius \(k\) and thickness \(\mathrm{d}k\). Such a shell has a volume \(\mathrm{d}V = \frac18 (4 \pi k^2) \mathrm{d}k\), where the factor \(\frac18\) is because we’re only considering an octant of the k-space. As each state occupies a volume \(\pi^3/V\) in k-space, and we have two electrons per unit volume, the number of states in our shell is given by
To get the total energy of the ground states, we multiply the number of states per shell with the energy per state, then integrate over all shells to get
The ground state energy \(E_\mathrm{tot}\) is the quantum-mechanical analog of the internal energy \(U\) of an (ideal) gas in (classical) thermodynamics. Unsurprisingly, it increases if we increase the number of electrons at constant volume, and decreases if we increase the volume at constant number of particles. Changing the volume of the solid will thus require work; therefore, the free electron gas must have an (internal) pressure. We can easily calculate this pressure using the known relations from thermodynamics. The first law of thermodynamics states that changes in the internal energy can only come from a flow of heat or by doing work on the gas, \(\Delta U = \Delta Q + \Delta W\). We are considering the ground state of our solid, i.e., the lowest possible energy state. In the language of thermodynamics, this state would be the coldest possible state, or \(T=0\). Work can be done by changing the volume, according to the relation \(\mathrm{d}W = - p \mathrm{d}V\), with \(p\) the pressure. We can then calculate the pressure from the relation
Equating the internal energy with the ground state energy, we obtain for the pressure
The relation \(p = \frac23 E_\mathrm{tot}/V\) is the same as the one we obtain for a classical ideal gas. The origin of the pressure however is completely different: in a classical ideal gas, the pressure on the walls of the container is due to the kinetic energy (and thus the temperature) of the gas particles. Here, we’ve explicitly set the temperature to zero, which classically would reduce the energy to zero as well. However, due to the combination of the impossibility of a zero-energy state in quantum mechanics and the Pauli exclusion principle which ‘forces’ electrons in higher energy states even in the lowest-energy configuration, we still get a nonzero pressure[12]. This internal pressure is therefore due to purely quantum-mechanical effects. This pressure counteracts the attractive forces that keep the solid together, leading to a stable state in which the solid becomes a ‘hard’ object, that neither expands nor collapses, and interacts (as a whole) with other solid objects. While compressing the solid is very difficult, solids do support longitudinal waves (aka sound), which are essentially waves of local compression and expansion. These waves carry energy and can be quantized; the corresponding quantum ‘particles’ are known as phonons.
5.3.2. Nuclear interactions and band structure#
The free electron gas gives us the Fermi surface and corresponding stabilizing internal pressure of a solid, but cannot distinguish between electrically conducting and isolating materials. To account for those, we need to add a layer of detail: interaction of the electrons with the nuclei. Fortunately, even a very crude model will do to generate the key concepts. Two popular choices include the Dirac comb, in which each nucleus is represented by a Dirac delta function potential, and the Kronig-Penney model, in which each nucleus is a finite rectangular potential. Between the nuclei, the potential remains zero. Qualitatively the results will be the same, as they both represent periodic potentials. Of course, a truly periodic potential needs to continue indefinitely, and we’ve seen in the previous section that the interactions of the electrons with the boundaries matter. Here however, we’re going to focus on the region between two nuclei, so in contrast to the global picture we studied above, we’re working with a local picture now. Compared to the typical spacing \(a\) between two nuclei, we can treat the solid as infinite; moreover, we will already get the qualitative results from a one-dimensional picture.
The key concept once again is going to be symmetry: if our potential is periodic, shifting everything by one (or multiple) periods shouldn’t change anything. The ‘anything’ is going to be ‘anything we can measure’, i.e., the square of the wave function.
Theorem 5.2 (Bloch’s theorem)
For a periodic potential with period \(a\), i.e., \(V(x) = V(x+a)\) for all \(x \in [0,a]\), the solutions to the time-independent Schrödinger equation satisfy
for some constant \(K\).
Proof. We define the displacement operator \(\hat{D}\) as \(\hat{D} f(x) = f(x+a)\). If the potential \(V(x)\) is periodic, then this operator commutes with the Hamiltonian, i.e. \([\hat{D}, \hat{H}] = 0\). Therefore, \(\hat{H}\) and \(\hat{D}\) share a complete common set of eigenfunctions. If \(\psi(x)\) is such an eigenfunction, it is a solution of the Schrödinger equation, and, as an eigenfunction of \(\hat{D}\), satisfies
For a normalizable function \(\psi(x)\), the eigenvalue \(\lambda\) cannot be zero, as that would immediately give \(\psi(x) = 0\). Therefore we can always write \(\lambda = \exp(i K a)\) for some \(K\), which gives equation (5.44).
Working in the local picture, we take our crystal lattice to be infinite, which we can do by imposing periodic boundary conditions on our wave function \(\psi(x)\):
where \(N\) is again the number of lattice sites (a very large number). By Bloch’s theorem we then have
so \(\exp(i N K a) = 1\), and \(NKa\) must be a multiple of \(2\pi\), or
Equation (5.45) makes our life much easier: rather than having to solve the Schrödinger equation for our entire material, we only need to solve it in a single ‘cell’ (an interval between two nuclei); we get every other cell for free by applying equation (5.44) with \(K\) given by (5.45). For the solution in the single cell, we need to choose a potential; as we already know how to deal with Dirac peaks, we’ll work with the Dirac comb here:
The general solution to the corresponding Schrödinger equation in the cell between the nucleus at \(x=0\) and the nucleus at \(x=a\) is the familiar combination of sines and cosines
where \(k = \sqrt{2mE}/\hbar\), as before. Using the periodicity of the potential, we can now also write down the solution in the cell between the nucleus at \(x=-a\) and the nucleus at \(x=0\):
We can find the values of the constants \(A\) and \(B\) by imposing boundary conditions at \(x=0\). Continuity of \(\psi(x)\) gives
Because the potential is infinite at \(x=0\), the derivative of \(\psi(x)\) makes a jump there (see Section 2.2.3):
Using equation (5.46) to express \(A\) in terms of \(B\) and substituting in (5.47), we get an equation in which each term is proportional to \(B\) (as it should be, because one of the integration constants should be set by normalization), and we are left with a single equation for the single unknown parameter \(k\) (and thus the energy):
We cannot solve equation (5.48) analytically. Graphically however we can easily see where the solutions lie. To that end, we define \(z=ka\) and \(\beta = m \alpha a / \hbar^2\), so we can re-write equation (5.48) as
The function \(f(z)\) is plotted in Fig. 5.9(a). The right-hand side of equation (5.49) is limited to the interval \([-1, 1]\); therefore, there can be no solutions if \(f(z)\) is outside this interval. Because \(N\) is very large and \(n\) can take any integer value, there are many solutions when \(f(z)\) lies between \(-1\) and \(1\). These solutions represent a band of energies that electrons in the solid can all occupy, whereas values of \(z\) for which \(|f(z)|>1\) represent a gap in which no solutions exist. The resulting electronic structure is known as the band-gap structure, see Fig. 5.9(b). The properties of the material now depend on how many states of the band are occupied if the electrons are in the ground state. If a band is partially filled, exciting an electron is very easy, as the energy levels in the band are close together: such a material would be an electrical conductor. On the other hand, if in the ground state all levels in the highest band are filled, exciting an electron would require a large amount of energy, as it would have to jump the gap to the next band. Such materials represent electrical insulators.
Fig. 5.9 Conductance bands and insulating gaps in solid materials. (a) Graphical solution of equation (5.49). Solutions only exist in the gray regions where \(|f(z)| \leq 1\). (b) Conductance bands alternated by gaps in the energy levels of a solid. If in the ground state all levels right up to a gap are filled, exciting an electron is difficult and the material behaves as an insulator. If a band is only partially filled, exciting an electron is easy and the material is a conductor.#
\newcounter{varprinciplesectioncounter} \setcounter{varprinciplesectioncounter}{\value{section}}
5.4. Appendix: Elliptical coordinates#
Elliptic and elliptical coordinates are two orthogonal coordinate systems (in two and three dimensions, respectively), which turn out to be immensely useful when calculating the integrals in this chapter. The most useful definition for our purposes[13] involves the distances \(r_1\) and \(r_2\) to two reference points, which we’ll put one the \(z\)-axis, one at \(-\frac12 R\) and one at \(+\frac12 R\) (these could be the nuclei in a diatomic molecule), see Fig. 5.10. This system is by construction symmetric to rotations about the \(z\)-axis, so let us introduce (for the full 3D case) a rotation angle \(\phi\) around the \(z\) axis, starting on the \(x\)-axis, just like in cylindrical or spherical coordinates. For our other two coordinates we define
where \(1 \leq \mu < \infty\) and \(-1 \leq \nu \leq 1\) (as you can easily check for yourself). It is not difficult to relate the \((\mu, \nu, \phi)\) coordinates to Cartesian and cylindrical ones. Consider a point \(\bm{r} = (x, y, z) = (\rho \cos\phi, \rho \sin\phi, z)\) in three-dimensional space (with \(\rho\) the distance of the point to the \(z\)-axis, and \(\phi\) the angle with the \(x\)-axis), then
from which we can solve for \(z\) and \(\rho\):
We have \(x = \rho \cos\phi\), \(y = \rho \sin\phi\) as in cylindrical coordinates. Inversely, we can also substitute the expressions for \(r_1\) and \(r_2\) in terms of the Cartesian coordinates into the (5.50) to get explicit expressions for \(\mu\) and \(\nu\) in terms of \(x\), \(y\), and \(z\); as for cylindrical coordinates, we also have \(\tan \phi = y/x\).
Fig. 5.10 Definition of our elliptical coordinates. Starting from two points at \(\pm \frac12 R\) on the \(z\) axis, we can calculate the distances \(r_1\) and \(r_2\) to both points for any point in space (red dot). Points for which the sum of \(r_1\) and \(r_2\) are equal define an ellipse (blue line) in the \((z, \rho)\) plane. Points for which the difference between \(r_1\) and \(r_2\) are equal define a hyperbola (orange line) in the \((z, \rho)\) plane. The ellipse and hyperbola have one intersection point for positive values of \(\rho\). By adding a rotation angle \(\phi\) around the \(z\)-axis, we can uniquely identify any point in three-dimensional space by specifying the ellipse (\(\mu\)), hyperbola (\(\nu\)) and angle (\(\phi\)).#
The coordinates defined by equation (5.50) are called elliptic (in the two-dimensional \((z, \rho)\) plane, or \((z, x)\) plane if we take \(\phi = 0\)) because the lines of constant \(\mu\) form ellipses, with foci at \(z = \pm \frac12 R\). To see that this is true, we re-write equation (5.51), by first taking the square of each equation, and then using the first in the second to eliminate \(\nu\), which gives:
which we can easily rewrite to
Equation (5.52) defines an ellipse for every value of \(\mu > 1\) (for \(\mu = 1\) we have the line from \(-\frac12 R\) to \(+\frac12 R\) on the \(z\)-axis). Likewise, we could have eliminated \(\mu\) instead of \(\nu\), which gives
so the lines of constant \(\nu\) for \(-1 \leq \nu \leq 1\) are hyperbolae with foci at \(z = \pm \frac12 R\).
Fig. 5.11 Lines of constant \(\mu\) (blue ellipses) and \(\nu\) (orange hyperbolae) in the \((z, \rho)\) plane.#
Figuring out the area and volume elements of our elliptic and elliptical coordinates is a straightforward exercise in coordinate transformations. For the area element, we find
For the volume element, we start with Cartesian coordinates and first make a substitution to cylindrical ones, which gives us an easy and surprisingly succinct result:
In these elliptical coordinates, the various integrals of the LCAO, especially for diatomic molecules, are much easier to evaluate than in spherical or cylindrical coordinates.
5.5. Appendix: Molecular integrals#
The following integrals appear in the calculation of the energies of diatomic molecules; the first three also appear in the expressions for the hydrogen molecule ion.
Overlap integral \(S\) (equation (5.14)):
Direct or Coulomb integral \(D\) (equation (5.16)):
Exchange integral \(X\) (equation (5.17)):
The repulsion between the charge density of electron 1 on nucleus A and the charge density of electron 2 on nucleus B, equation (5.33)a:
The repulsion between the overlap charge density of electron 1 and the overlap charge density of electron 2, equation (5.33)b:
where
with \(s=R/a\), \(S(s)\) the overlap integral given by equation (5.56), \(\gamma \; (= 0.57722...)\) is Euler’s constant, and \(\mathrm{Ei}(x)\) the exponential integral defined as
The repulsion between the charge density of electron 1 on nucleus A and and the overlap charge density of electron 2, equation (5.33)c:
The repulsion between the charge density of electron 1 on nucleus A and and the charge density of electron 2, also on A, equation (5.33)d:
5.6. Problems#
Exercise 5.1 (Estimating the electron interaction energy in the ground state of helium)
In Section 5.1.1, we discussed that while we can solve for the eigenfunctions of the electron in the hydrogen atom, we cannot do so for the two electrons in the helium atom, due to their mutual repulsion. As a zeroth-order approximation, we can just ignore this interaction, and the ground state of helium (with two electrons at \(\bm{r}_1\) and \(\bm{r}_2\) simply becomes the product of two ground states of hydrogen, as given by equation (5.5).
Show that the energy of the helium ground state for two non-interacting electrons, as approximated by equation (5.5), would be \(8\) times the ground state energy of hydrogen. (The corresponding value is \(-109\;\mathrm{eV}\), while the actual (experimental) value is \(-79\;\mathrm{eV}\), so the approximation is really quite bad). We now wish to calculate \(\Braket{(1/|\bm{r}_1 - \bm{r}_2|)}\) for the product state of equation (5.5). There are (at least) two ways to do this. One is to use spherical coordinates for both \(\bm{r}_1\) and \(\bm{r}_2\); you can then set the polar axis of \(\bm{r}_2\) along \(\bm{r}_1\), so that \(|\bm{r}_1 - \bm{r}_2| = \sqrt{r_1^2 + r_2^2 - 2 r_1 r_2 \cos\theta}\). Although this is doable, it gives you complicated integrals. A ‘cleaner’ method is to make use of the symmetry of the problem. We’re integrating over all space twice, which we can interpret as integrating over the positions of the two electrons, but also as integrating over the separation vector \(\bm{R} = \bm{r}_1 - \bm{r}_2\) between the electrons and the position \(\bm{r}\) of the nucleus. We can then use elliptical coordinates to re-write the distances \(r_1\) and \(r_2\) between the endpoints of \(\bm{R}\) and the position \(\bm{r}\) of the nucleus in an easily integrable form, see Section 5.4.
Write the integral for \(\Braket{(1/|\bm{r}_1 - \bm{r}_2|)}\) in terms of \(\bm{R}\) and \(\bm{r}\). Convert to spherical coordinates for the \(\bm{R}\) integral.
Write the integral for \(\bm{r}\) in elliptical coordinates. You should have three non-trivial integrals left (over \(R\), \(\mu\), and \(\nu\)).
Evaluate each of the integrals in turn. You may use that
\[\int_0^\infty x^n e^{-x/a} \mathrm{d}x = n! \, a^{n+1},\]for any nonnegative integer \(n\). The final answer (which you’ll need for (e)) is \(\frac{5}{4a}\).
Use the result in (d) to estimate the electron interaction energy in the ground state of helium. Express your answer in electron volts, and add it to the energy of (a) to get a corrected estimate of the ground-state energy. Compare your answer to the experimental value.
Exercise 5.2 (The LCAO for the molecular hydrogen ion)
Derive the expressions for the direct (\(D\)) and exchange (\(X\)) integrals given in equations (5.16) and (5.17).
Derive the expression (5.20) for the expectation value of the molecular hydrogen ion Hamiltonian (equation (5.11)) in the antisymmetric LCAO state (equation (5.19)).
Re-derive the equations for the energy of the symmetric and antisymmetric LCAO of the hydrogen ion from the general condition (5.27) in Section 5.2.3. Hint: first express the matrix elements \(H_{ik}\) and \(S_{ik}\) in terms of the integrals \(D\), \(X\) and \(S\) of Section 5.2.2; you’ll get two symmetric \(2 \times 2\) matrices, after which solving for the energies from equation (5.27) will be easy.
Exercise 5.3 (Oscillatoions in the molecular hydrogen ion)
The functional form of the difference between the energy estimate of the molecular hydrogen ion and the hydrogen atom, as given by equation (5.18) and plotted in Fig. 5.2 (blue line) can be written as
where \(s = R/a\) and \(F\) is measured in units of \(-E_1\). As we’ve seen, this function has a minimum that is lower than zero, suggesting that covalent bonds can form. However, as we’ve also seen previously, particles (here the nuclei) in a potential well are never stationary, and their ground state energies are larger than the minimum. In this problem we’ll estimate those energies, and show that the nuclei remain safely bound.
Numerically find the position and value of the minimum of \(F(s)\).
As we can locally approximate \(F(s)\) about its minimum with a parabola, we can estimate the frequency \(\omega\) of the corresponding harmonic oscillator potential from the second derivative of the function \(F(s)\). (If you don’t remember how to do so from classical mechanics, consider the comparison with a classical harmonic oscillator with potential \(V(x) = \frac12 k x^2\)).
Show that the ground state energy of the harmonic oscillator potential approximation (i.e., \(\frac12 \hbar \omega\)) does not exceed the binding energy of the two protons.
Estimate how many bound vibrational modes the molecular hydrogen ion has.