A schematic representation of the ingredients of a molecule: nuclei (black) and electrons (red).

Quantum theory of molecular many-particle systems

Post by: Bjoern Baumeier
December 24, 2018
No Comment

I. Quantum theory of molecular many-particle systems

These notes are part of a full lecture about “Charge and Energy Transfer Dynamics in Molecular Systems“. See its full Table of Contents.

Molecular Schrödinger Equation
Born-Oppenheimer Separation and the adiabatic approximation
Electronic structure methods: overview

Non-interacting single-particles
The Hartree-Fock method
Density-Functional Theory
Comparison of single-particle and excitation energies

Fermi’s Golden Rule
Harmonic approximation for nuclear motion

I.1 Molecular Schrödinger Equation

For the moment, let us consider a rather general picture of a molecule as indicated in Figure ??: an aggregate structure consisting of $N$ nuclei and $n$ electrons. The coordinates $\mathbf{R}_\alpha$ (spins $\Sigma_\alpha\mathrm{)}$ of the individual nuclei and $\mathbf{r}_i$ $\mathrm{(}\sigma_i\mathrm{)}$ of the individual electrons are combined into the supervariables $\mathbf{R}=\left(\mathbf{R}_1, \ldots , \mathbf{R}_N\right)$ $\mathrm{(}\Sigma = \left(\Sigma_1, \ldots , \Sigma_N\right) \mathrm{)}$ and $\mathbf{r}=\left(\mathbf{r}_1, \ldots , \mathbf{r}_n\right)$ $\mathrm{(}\sigma = \left(\sigma_1, \ldots , \sigma_n\right) \mathrm{),}$ respectively.

The full quantum-mechanical information about the molecule is — at least in principle — obtained by solving the time-dependent Schrödinger equation

(1) $\begin{equation*} \hat{H}_\mathrm{mol} \Psi(\mathbf{r}\sigma, \mathbf{R}\Sigma,t) = i \hbar \frac{\partial}{\partial t} \Psi(\mathbf{r}\sigma, \mathbf{R}\Sigma,t) \end{equation*}$

The Hamiltonian of the molecule is (in the absence of external fields) given by

(2) $\begin{eqnarray*} \hat{H}_\mathrm{mol} &=& \hat{T}_\mathrm{nuc} + \hat{V}_\mathrm{nuc-nuc} + \hat{T}_\mathrm{el} + \hat{V}_\mathrm{el-el} + \hat{V}_\mathrm{nuc-el}\\ &=&\sum_{\alpha=1}^N{\frac{\hat{P}_\alpha^2}{2M_\alpha}} + \frac{1}{2} \sum_{\alpha , \beta}^{N,N}{\frac{Z_\alpha Z_\beta e^2}{\vert \mathbf{R}_\alpha - \mathbf{R}_\beta \vert}}\\ &&+\sum_{i=1}^n{\frac{\hat{p}_i^2}{2m}} + \frac{1}{2} \sum_{i , j}^{n,n}{\frac{e^2}{\vert \mathbf{r}_i - \mathbf{r}_j \vert}} \\ &&-\sum_{\alpha=1}^N{\sum_{i=1}^n{\frac{Z_\alpha e^2}{\vert \mathbf{r}_i-\mathbf{R}_\alpha \vert}}} \end{eqnarray*}$

The Hamiltonian in Eq.2

is except for relativistic corrections exact;
does not explicitly depend on the time $t\mathrm{;}$
is exactly solvable only for $N=n=1$ (hydrogen atom)
“contains” all phenomena, such as
- structural properties (equilibrium geometries)
- dynamical properties (vibrations)
- electronic properties (ionization spectra, UPS)
- optical properties (absorption, emission)
- transport properties
- …

From ii. it follows that wave function $\Psi(\mathbf{r}\sigma, \mathbf{R}\Sigma,t)$ as the solution of the time-dependent Schrödinger equation Eq.1 factorizes as

(3) $\begin{equation*} \Psi(\mathbf{r}\sigma, \mathbf{R}\Sigma,t) = \psi(\mathbf{r}\sigma, \mathbf{R}\Sigma)\cdot u(t) \end{equation*}$

where $u(t)=\exp{\left(-i\frac{E_\mathrm{mol}}{\hbar}t\right)}$ and $\psi(\mathbf{r}\sigma, \mathbf{R}\Sigma)$ and $E_\mathrm{mol}$ are obtained as solutions of the stationary Schrödinger equation

(4) $\begin{equation*} \hat{H}_\mathrm{mol} \psi(\mathbf{r}\sigma, \mathbf{R}\Sigma) =E_\mathrm{mol} \psi(\mathbf{r}\sigma, \mathbf{R}\Sigma) \end{equation*}$

The task that we face in describing charge and energy transfer processes is clear: solve Eq.4 and extract the relevant properties from the solution. Obviously, this is impossible to achieve exactly and instead has to rely on physically sensible approximations.

Back to Table of Contents.

I.2 Born-Oppenheimer Separation and the adiabatic approximation

The first approximation to introduce deals with the relevant time scales of nuclear and electronic motion. Two rather simple arguments can be put forward to motivate this:

We know that stable molecules exist!
From this very basic realization, it follows that the forces acting on the electrons and the nuclei must be similar in this stable equilibrium:
(5) $\begin{eqnarray*} \mathbf{F}_i = m \frac{\partial^2 \mathbf{r}_i}{\partial t^2}= M_\alpha \frac{\partial^2 \mathbf{R}_\alpha}{\partial t^2} = \mathbf{F}_\alpha\\ \Rightarrow \frac{\partial^2 \mathbf{R}_\alpha}{\partial t^2} \approx \frac{m}{M_\alpha}\frac{\partial^2 \mathbf{r}_i}{\partial t^2} \approx 10^{-4} \frac{\partial^2 \mathbf{r}_i}{\partial t^2} \end{eqnarray*}$

As a consequence, the dynamics of the nuclei is much slower than that of the electrons. In other words, the electrons adjust instantaneously to the nuclear motion, i.e., the electrons move adiabatically.

To express this situation in formal terms, we consider a fixed arrangement of nuclei $\mathbf{R}\mathrm{.}$ The electronic Hamiltonian representing the electronic system that interacts with the fixed nuclear configuration reads

(6) $\begin{equation*} \hat{H}_\mathrm{el} = \hat{H}_\mathrm{el}(\mathbf{R}) = \hat{T}_\mathrm{el} + \hat{V}_\mathrm{el-el} + \hat{V}_\mathrm{nuc-el}(\mathbf{R}). \end{equation*}$

In this situation, $\mathbf{R}$ is no longer a variable of the electronic system, but a fixed parameter for the electronic degrees of freedom. The corresponding stationary electronic Schrödinger equation is given by

(7) $\begin{equation*} \hat{H}_\mathrm{el}(\mathbf{R}) \Phi_\nu(\mathbf{r};\mathbf{R}) = E_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R}), \end{equation*}$

where $\left\{\Phi_\nu(\mathbf{r};\mathbf{R})\right\}$ is a set of adiabatic electronic wave functions. Those can be used as a basis to expand the molecular wave function $\Psi(\mathbf{r},\mathbf{R})$ according to

(8) $\begin{equation*} \psi(\mathbf{r},\mathbf{R}) = \sum_\nu{\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})} . \end{equation*}$

Entering this Born-Oppenheimer separated wave function into Eq.4 yield

(9) $\begin{eqnarray*} \hat{H}_\mathrm{mol} \psi(\mathbf{r}, \mathbf{R}) &=& \left( \hat{H}_\mathrm{el}(\mathbf{R}) + \hat{T}_\mathrm{nuc} + \hat{V}_\mathrm{nuc-nuc} \right) \sum_\nu{\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})} \\ &=& E_\mathrm{mol} \psi(\mathbf{r}\sigma, \mathbf{R}\Sigma) \\ &=& \sum_\nu{ \left( E_\nu(\mathbf{R}) + \hat{T}_\mathrm{nuc} + \hat{V}_\mathrm{nuc-nuc} \right) \chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})} \end{eqnarray*}$

If we now first multiply this equation by $\Phi_\mu^*(\mathbf{r};\mathbf{R})$ and then integrate over the electronic coordinates, we obtain an expression for the nuclear functions $\chi_\nu(\mathbf{R})$

(10) $\begin{equation*} E_\mathrm{mol}\chi_\nu(\mathbf{R}) = \left(E_\nu(\mathbf{R}) + \hat{V}_\mathrm{nuc-nuc} \right)\chi_\nu(\mathbf{R}) +\int{ \Phi_\mu^*(\mathbf{r};\mathbf{R}) \hat{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})\,d\mathbf{r}}. \end{equation*}$

One has to pay special attention evaluating the action of the kinetic energy operator of the nuclei $\hat{T}_\mathrm{nuc}$ on the terms of the adiabatic wave function, marked red, within the integral above. Using the product rule, one finds

(11) $\begin{eqnarray*} \hat{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R}) &=& \sum_{\alpha=1}^N{-\frac{\hbar^2}{2M_\alpha}\nabla_{\mathbf{R}_\alpha}^2 \chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R}) }\\ &=&\Phi_\nu(\mathbf{r};\mathbf{R}) \hat{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\\ &&\underbrace{-2\sum_{\alpha=1}^N{\frac{\hbar^2}{2M_\alpha}\frac{\partial}{\partial \mathbf{R}_\alpha}\Phi_\nu(\mathbf{r};\mathbf{R}) \frac{\partial}{\partial \mathbf{R}_\alpha} \chi_\nu(\mathbf{R}) } + \chi_\nu(\mathbf{R})\hat{T}_\mathrm{nuc}\Phi_\nu(\mathbf{r};\mathbf{R})}_{\bar{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})} \\&=&\Phi_\nu(\mathbf{r};\mathbf{R}) \hat{T}_\mathrm{nuc}\chi_\nu(\mathbf{R}) + \bar{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R}) \end{eqnarray*}$

We can now enter this result into the integral in Eq.10

(12) $\begin{equation*} \int{ \Phi_\mu^*(\mathbf{r};\mathbf{R}) \hat{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})\,d\mathbf{r}} = \hat{T}_\mathrm{nuc}\chi_\nu(\mathbf{R}) + \int{ \Phi_\mu^*(\mathbf{r};\mathbf{R}) \bar{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})\,d\mathbf{r}} \end{equation*}$

The remaining integral

(13) $\begin{equation*} A_{\nu\mu} = \int{ \Phi_\mu^*(\mathbf{r};\mathbf{R}) \bar{T}_\mathrm{nuc}\chi_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R})\,d\mathbf{r}} \end{equation*}$

defines matrix elements of the transition between electronic states $\nu$ and $\mu$ induced by the dynamics of the nuclei. This means that in order to determine the functions $\chi_\nu(\mathbf{R})$ one has to solve the coupled set of equations

(14) $\begin{equation*} E_\mathrm{mol}\chi_\nu(\mathbf{R}) = \left(E_\nu(\mathbf{R}) +\hat{T}_\mathrm{nuc} + \hat{V}_\mathrm{nuc-nuc} \right)\chi_\nu(\mathbf{R}) + \sum_\mu{A_{\nu\mu}\chi_\mu(\mathbf{R})} \end{equation*}$

In the adiabatic approximation it is assumed that $A_{\mu\nu}=0\mathrm{,}$ i.e., there are no transitions between different electronic states, and the nuclear motion for each electronic state $\nu$ is determined by

(15) $\begin{equation*} E_\mathrm{mol}\chi_\nu(\mathbf{R}) = \left(E_\nu(\mathbf{R}) +\hat{T}_\mathrm{nuc} + \hat{V}_\mathrm{nuc-nuc} \right)\chi_\nu(\mathbf{R}) . \end{equation*}$

This describes the motion of the nuclei in an effective potential

(16) $\begin{equation*} U_\nu(\mathbf{R}) = E_\nu(\mathbf{R}) + \hat{V}_\mathrm{nuc-nuc}(\mathbf{R}) , \end{equation*}$

which is also called the potential energy surface (PES). Note that for each electronic states $\nu\mathrm{,}$ there is a unique PES associated with it, as indicated in Fig.??.

1D representation of a potential energy surfaces for different electronic states of a molecule: neutral (black), optically excited (blue), and cation (red).

In practical calculations, one has to

Fix nuclear coordinates $\mathbf{R}$ and solve the electronic problem in Eq.7, obtaining $E_\nu(\mathbf{R})$ and $\Phi_\nu(\mathbf{r};\mathbf{R})\mathrm{.}$
Solve Eq.15 with the effective potential $U_\nu(\mathbf{R})\mathrm{.}$

Back to Table of Contents.

I.3 Electronic structure methods: an overview

In this part, we will take a look at fundamental approaches for the solution of the electronic problem commonly used in quantum chemistry. It is obvious that at this point we cannot give a fully detailed account of those methods and refer to standard textbooks (e.g., C.J. Kramer, “Essentials of Computational Chemistry”, 2nd Edition, Wiley, 2008) for an in-depth discussion. The purpose of this part is instead to convey the basic ideas, introduce typically used terminology, and illustrate the strength and weaknesses for the various approaches.

What we seek are the solutions of the $n$ electron Schrödinger equation 7

(17) $\begin{equation*} \hat{H}_\mathrm{el}(\mathbf{R}) \Phi_\nu(\mathbf{r};\mathbf{R}) = E_\nu(\mathbf{R})\Phi_\nu(\mathbf{r};\mathbf{R}), \end{equation*}$

where the $\Phi_\nu(\mathbf{r};\mathbf{R})$ are adiabatic wave functions of $n$ interacting electrons and the corresponding Hamiltonian is

(18) $\begin{eqnarray*} \hat{H}_\mathrm{el}(\mathbf{R}) &=& \hat{T}_\mathrm{el} + \hat{V}_\mathrm{el-el} + \hat{V}_\mathrm{nuc-el}(\mathbf{R})\\ &=& \sum_{i=1}^n{\left(\frac{\hat{p}_i^2}{2m}-\sum_{\alpha=1}^N{\frac{e^2Z_\alpha}{\vert \mathbf{r}_i-\mathbf{R}_\alpha\vert}}\right)} +\frac{1}{2}\sum_{i,j}^{n,n}{\frac{e^2}{\vert \mathbf{r}_i - \mathbf{r}_j\vert}}\\ &=&\sum_{i=1}^n{\hat{h}_i} + \sum_{i,j}^{n,n}{v(\mathbf{r}_i,\mathbf{r}_j)} . \end{eqnarray*}$

One can clearly see that the electronic Hamiltonian comprises a sum of single-particle Hamiltonians $\hat{h}_i$ and a many-body Coulomb term $v\mathrm{.}$ The various approaches presented in the following essentially all differ in the approximations made to take this many-body interaction into account.

Non-interacting single-particles

The simplest possible approximation that one can make is the limit of vanishing Coulomb interactions, or non-interacting single-particles. In this case, the electronic Hamiltonian consists solely of the sum of single-particle terms

(19) $\begin{equation*} \hat{H}_\mathrm{el}(\mathbf{R}) = \sum_{i=1}^n{\hat{h}_i} \end{equation*}$

and the corresponding $n$ electron wave function $\Phi_\nu^0(\mathbf{r};\mathbf{R})$ is simply a product of single-particle functions

(20) $\begin{equation*} \Phi_\nu^0(\mathbf{r};\mathbf{R}) = \prod_{i=1}^n{\phi_{\nu_i}^0(\mathbf{r}_i)} . \end{equation*}$

These single-particle functions are solutions to the single-particle Schrödinger equation

(21) $\begin{equation*} \hat{h}_i \phi_{\nu_i}^0(\mathbf{r}_i) = \varepsilon_{\nu_i}^0 \phi_{\nu_i}^0(\mathbf{r}_i) \end{equation*}$

where $\varepsilon_{\nu_i}^0$ is a single-particle energy.

What we have now managed is to transform a single equation for $n$ (non-interacting) particles to $n$ equations of single-particles, which are obviously much easier to treat. The total energy of the non-interacting particle system is given by

(22) $\begin{equation*} E_\nu^0=\sum_{i=1}^n{\epsilon_{\nu_i}^0}, \end{equation*}$

i.e., a sum of the respective single-particle energies. Note that $\nu$ denotes a set of single-particle quantum-numbers $\{\nu_i\}\mathrm{.}$ When forming the product function in Eq.20 one only has to take the Pauli exclusion-principle into account.

Filling up the single-particle states with electrons starting from the lowest energiesobserving the Pauli exclusion-principle. — Filling up the single-particle states with electrons starting from the lowest energies observing the Pauli exclusion-principle.

One can easily form the configuration of minimal $E_\nu^0=E_0^0$ , i.e. the ground state energy by sorting all $\epsilon_{\nu_i}$ according to their value and the filling up each energy level starting from the bottom with two electrons of opposite spin, see Fig.??. The single-particle function $\phi^0(\mathbf{r}_i)$ of the last filled energy level is called the highest occupied molecular orbital (HOMO), the one of the first empty level the lowest unoccupied molecular orbital (LUMO).

The beauty of this approach is that one can straightforwardly relate important quantities from the non-interacting single-particle energies:

1. electron removal
  If we remove the $l\mathrm{-th}$ electron, the energy of the $n-1_l$ electron system is
  
  (23) $\begin{eqnarray*} E^0(n-1_l)=\sum_{i=1}^n{\epsilon_{i}^0}-\epsilon_{l} = E^0(n)-\epsilon_{l}\\ \Rightarrow E^0(n-1_l)-E^0(n)=-\epsilon_{l}=\mathrm{IP}_l \end{eqnarray*}$
  
  The occupied single-particle energies are therefore the negative of the respective ionization energy.
2. electron addition
  If we add an electron to the $k\mathrm{-th}$ energy level, the energy of the $n+1_k$ electron system is
  
  (24) $\begin{eqnarray*} E^0(n+1_k)=\sum_{i=1}^n{\epsilon_{i}^0}+\epsilon_{k} = E^0(n)+\epsilon_{k}\\ \Rightarrow E^0(n+1_k)-E^0(n)=\epsilon_{k}=-\mathrm{EA}_k \end{eqnarray*}$
  
  The unoccupied single-particle energies are therefore the negative of the respective electron affinities (though the definition of the sign may differ).
3. electron promotion

If we promote an electron from an occupied level $l$ to the previously empty $k\mathrm{-th}$ energy level (e.g. by optical excitation), the total energy of is

(25) $\begin{equation*} E^0(n-1_l+1_k)= E^0(n)+\epsilon_{k}-\epsilon_{l} \end{equation*}$

So the, e.g. energy of absorption, is equal to the difference in single-particle energies.

Often, basic details of the electronic structure of molecules related to charge and energy transfer processes are discussed in term of this single-particle picture. However, the downside is that, of course, the electrons of realistic molecular systems do interact. Still, one can certain within limits stay within this framework taking the interactions of the electrons into account by transforming the full interacting problem into what is called an effective single-particle problem, as will be shown in the following.

The Hartree-Fock method

One aspect that we have not addressed in the previous section is the explicit Fermion nature of the electrons that requires the $n$ electron wave function $\Phi_\nu(\mathbf{r};\mathbf{R})$ to be antisymmetric with regards to particle exchange. This requirement is not fulfilled by the simple product ansatz in Eq.20. One can, however, construct an appropriate function from single-particle functions using a determinant ansatz

(26) $\begin{equation*} \Phi_\nu(\mathbf{r};\mathbf{R}) = \frac{1}{\sqrt{n!}} \left\vert \begin{array}{ccc} \phi_1(\mathbf{r}_1) & \ldots & \phi_1(\mathbf{r}_n)\\ \vdots & \ddots & \vdots \\ \phi_n(\mathbf{r}_1) & \ldots & \phi_n(\mathbf{r}_n) \end{array}\right\vert \end{equation*}$

often referred to as Slater determinant. Now, instead of starting from predetermined single-particle functions and enforcing antisymmetry, we can start from the requirement of antisymmetry and use the variational principle to derive a set of equations that determine suitable effective single-particles for the interacting case. This is the principle idea of the Hartree-Fock method. Let us focus on the ground state energy $E_O^\mathrm{HF}$ and consider a Slater determinant as in Eq.26. From the variational principle it holds that the energy as a functional of that determinant approximates the true ground state energy from above, i.e.,

(27) $\begin{equation*} E_0^\mathrm{HF} = E[\Phi^\mathrm{HF}] = \frac{\langle \Phi^\mathrm{HF} \vert \hat{H}_\mathrm{el} \vert\Phi^\mathrm{HF} \rangle}{\langle \Phi^\mathrm{HF}\vert \Phi^\mathrm{HF} \rangle} \ge E_0 . \end{equation*}$

The task is now to minimize this functional by means of a non-linear variation with respect to the (undetermined) single-particle functions $\phi_i^\mathrm{HF}$ under the condition that those functions are normalized (introduction of Langragian multipliers $\epsilon_i^\mathrm{HF}$ )

(28) $\begin{equation*} \delta \left\{\langle \Phi^\mathrm{HF} \vert \hat{H}_\mathrm{el} \vert\Phi^\mathrm{HF} \rangle - \sum_{i=1}^n{\epsilon_i^\mathrm{HF}\left(\langle \Phi^\mathrm{HF}\vert \Phi^\mathrm{HF} \rangle -1 \right)} \right\} =0 \end{equation*}$

This variation can be performed using a functional derivative

(29) $\begin{equation*} \frac{\delta}{\delta \phi_j^{*}(x)} \left\{E[\Phi^\mathrm{HF}] - \sum_{i=1}^n{\epsilon_i^\mathrm{HF} \left( \int{\phi_i^*(x)\phi_i(x)\,dx}-1}\right)} =0 \right\} \end{equation*}$

where we introduced $x$ as a variable containing $\mathbf{r}$ and $\sigma$ and suppressed the superscript $\mathrm{HF}$ of the single-particle functions. After some algebra (see textbooks), one arrives at a set of equations that allow to determine the $\phi_j(x)$

(30) $\begin{eqnarray*} \left\{ -\frac{\hbar^2}{2m}\Delta - \sum_{\alpha=1}^N{\frac{e^2 Z_\alpha}{\vert \mathbf{r}-\mathbf{R}_\alpha \vert}}\right\} \phi_j(x) + \left\{ \sum_{i=1}^n{\int{\phi_i^*(x')v(x,x')\phi_i(x')\,dx'}}\right\} \phi_j(x) \\ -\sum_{i=1}^n{\int{\phi_i^*(x')v(x,x')\phi_j(x')\phi_i(x)\,dx'}} = \epsilon_j^\mathrm{HF} \phi_j(x) \end{eqnarray*}$

Introducing the densities

(31) $\begin{eqnarray*} \rho(x')&=&\sum_{i=1}^n{\phi_i^*(x')\phi_i(x')}\\ \rho(x,x')&=& \sum_{i=1}^n{\phi_i^*(x')\phi_i(x)} \end{eqnarray*}$

the Hartree-Fock equations read

(32) $\begin{eqnarray*} \left\{ -\frac{\hbar^2}{2m}\Delta + V_\mathrm{ext}(x) + \int{\rho(x')v(x,x')\,dx'} \right\} \phi_j^\mathrm{HF}(x)\\ -\int{\rho(x',x)v(x',x)\phi_j^\mathrm{HF}(x')\,dx'} = \epsilon_j^\mathrm{HF} \phi_j^\mathrm{HF}(x) . \end{eqnarray*}$

While the first integral in Eq.32 corresponds to the classical Hartree integral of the Coulomb interaction

(33) $\begin{equation*} V_H(x)=\int{\rho(x')v(x,x')\,dx'} \end{equation*}$

the second integral defines the exchange potential in operator form

(34) $\begin{equation*} V_x(x)=\int{\rho(x',x)v(x,x')\ldots\,dx'} \end{equation*}$

that is a direct result of the quantum-mechanical antisymmetry condition of the wave function. The $n$ electron problem has thus been mapped on a set of effective single-particle problems with the Hartree-Fock potential

(35) $\begin{equation*} \hat{V}_\mathrm{HF} = V_\textrm{ext}(x) + V_H(x) + \hat{V}_x(x) . \end{equation*}$

The total energy of the ground state is, however, not simply the sum of the “single-particle energies” $\epsilon_j^\mathrm{HF}$ (remember that these were introduced as Lagrangian multipliers during the variation) due to a double counting of $i-j$ interactions in $\epsilon_i^\mathrm{HF}$ and $\epsilon_j^\mathrm{HF}\mathrm{.}$ Therefore the total energy reads

(36) $\begin{equation*} E_0^\mathrm{HF} = \sum_{i=1}^n{\epsilon_i^\mathrm{HF}} - \frac{1}{2} \left( E_H + E_x \right) \end{equation*}$

where

(37) $\begin{eqnarray*} E_H &=& \sum_{i,j}^{n,n}{\int{ \phi_i^*(x')\phi_j^*(x'')v(x',x'') \phi_i(x')\phi_j(x'')\,dx'\,dx''}} \\ E_x &=& - \sum_{i,j}^{n,n}{\int{ \phi_i^*(x')\phi_j^*(x'')v(x',x'') \phi_j(x')\phi_i(x'')\,dx'\,dx''}} . \end{eqnarray*}$

In terms of the meaning of the $\epsilon_j^\mathrm{HF}$ , one finds

1. electron removal
  If we remove the $l\mathrm{-th}$ electron, the energy of the $n-1_l$ electron system is
  
  (38) $\begin{eqnarray*} E^\mathrm{HF}(n-1_l)=\sum_{i=1}^n{\epsilon_i^\mathrm{HF}}-\epsilon_{l}^\mathrm{HF}= E^\mathrm{HF}(n)-\epsilon_{l}^\mathrm{HF}\\ \Rightarrow E^\mathrm{HF}(n-1_l)-E^\mathrm{HF}(n)=-\epsilon_{l}^\mathrm{HF}=\mathrm{IP}_l \end{eqnarray*}$
  
  The occupied single-particle energies are therefore the negative of the respective ionization energy. This is Koopman’s theorem.
2. electron addition
  If we add an electron to the $k\mathrm{-th}$ energy level, the energy of the $n+1_k$ electron system is
  
  (39) $\begin{eqnarray*} E^\mathrm{HF}(n+1_k)=\sum_{i=1}^n{\epsilon_{i}^\mathrm{HF}}+\epsilon_{k} ^\mathrm{HF}= E^\mathrm{HF}(n)+\epsilon_{k}^\mathrm{HF}\\ \Rightarrow E^\mathrm{HF}(n+1_k)-E^\mathrm{HF}(n)=\epsilon_{k}^\mathrm{HF}=-\mathrm{EA}_k \end{eqnarray*}$
  
  The unoccupied single-particle energies are therefore the negative of the respective electron affinities (though the definition of the sign may differ). However, since the unoccupied states do not contribute to the exchange term in the Hartree-Fock equation, the respective energies contain a spurious repulsive electronic self-interaction resulting from the Hartree term. For the occupied states, this is compensated by the self-interaction in the exchange term!
3. electron promotion

If we promote an electron from an occupied level $l$ to the previously empty $k\mathrm{-th}$ energy level (e.g. by optical excitation), one can show that the total energy is

(40) $\begin{equation*} E^\mathrm{HF}(n-1_l+1_k)= E^\mathrm{HF}(n)+\epsilon_{k}^\mathrm{HF}-\epsilon_{l}^\mathrm{HF} + \mathrm{O}(1/n)\end{equation*}$

So the, e.g. energy of absorption, can be approximated by the difference in single-particle energies (with the limitations mentioned above).

Density-Functional Theory

Instead of dealing with the many-body wave function $\vert\Phi_0\rangle$ directly, two theorems by Hohenberg and Kohn relate the ground state to the electron density:

The ground state $\phi_0$ is a one-to-one functional of the particle density $\rho(\mathbf{r})\mathrm{.}$ Here, $\mathbf{r}$ is only a single spatial coordinate!
The energy functional
(41) $\begin{eqnarray*} E[\rho] &=& \int{V_\mathrm{ext}(\mathbf{r})\rho(\mathbf{r})\,d^3r} + \left\langle \Phi \vert \hat{T}+\hat{V}_\mathrm{el-el}\vert\Phi\right\rangle\\ &=& \int{V_\mathrm{ext}(\mathbf{r})\rho(\mathbf{r})\,d^3r} + F[\rho(\mathbf{r})] \end{eqnarray*}$

obeys a variational principle with respect to the particle density $\rho(\mathbf{r})$ and is minimal for the ground state density $\rho_0(\mathbf{r})$ :

(42) $\begin{equation*} E_0 = E[\rho_0] \leq E[\rho] . \end{equation*}$

Again, we introduce single-particle wave functions $\phi_i(\mathbf{r})$ to express the density

(43) $\begin{equation*} \rho(\mathbf{r}) = \sum_{i=1}^n{\vert \phi_i(\mathbf{r})\vert^2} \end{equation*}$

and proceed by separating the classical Coulomb effects (the Hartree energy) from any quantum-mechanical exchange-correlation interactions of the electrons

(44) $\begin{equation*} E[\rho] = T[\rho] + \int{V_\mathrm{ext}(\mathbf{r})\rho(\mathbf{r})\,d^3r} + \int{\int{\frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{\vert \mathbf{r} - \mathbf{r}'\vert}\,d^3r}d^3r'} + E_\mathrm{xc}[\rho] . \end{equation*}$

Variation of $E[\rho]$ with respect to $\rho$ can be rewritten as a functional derivative with respect to $\phi_i^*(\mathbf{r})$ under the constraint of normalization

(45) $\begin{equation*} \int{\phi_i^*(\mathbf{r})\phi_i(\mathbf{r})\,d^3r} = 1 \end{equation*}$

(formally identical procedure as in Hartree-Fock) and we again obtain a set of effective single-particle equations, the Kohn-Sham equations

(46) $\begin{equation*} \left\{-\frac{\hbar^2}{2m}\Delta + V_\mathrm{ext}(\mathbf{r}) + V_H(\mathbf{r}) + \frac{\delta E_\mathrm{xc}}{\delta \rho} \right\} \phi_i^\mathrm{KS}(\mathbf{r}) = \epsilon_i^\mathrm{KS} \phi_i^\mathrm{KS}(\mathbf{r}) . \end{equation*}$

The practical problem of applying these equations to obtain a formally exact single-particle representation of the $n$ electron ground state lies in the fact that the exchange-correlation functional $E_\mathrm{xc}$ that contains all quantum interactions is unknown, ans with it also the exchange-correlation potential

(47) $\begin{equation*} V_\mathrm{xc}[\rho] = \frac{\delta E_\mathrm{xc}}{\delta \rho} \end{equation*}$

in Eq.46. In practice, one therefore has to resort to physically motivated approximations to $V_\mathrm{xc}$ , such as

Local-density approximation (LDA)
Here one assumes that the real electron density at a certain point $\mathbf{r}$ can be locally approximated by that of the free-electron gas, for which $V_\mathrm{xc}$ is known and can be parametrized.
Generalized-gradient approximation (GGA, PBE, …)
As an extension to the LDA, the local gradient of the electron density is taken additionally into account.
hybrid functionals (B3LYP, PBE0, …)
Here, a mix of local or semi-local functionals with exact exchange contributions (see Hartree-Fock) is used.

Irrespective of that actual functional, we can determine the ground state energy again by removing the electron-electron interaction double counting

(48) $\begin{equation*} E^\mathrm{DFT}[\rho] = \sum_{i=1}^n{\epsilon_i^\mathrm{KS}}- \int{\int{\frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{\vert \mathbf{r}-\mathbf{r}'\vert}d^3r}d^3r'} - \int{V_\mathrm{xc}[\rho]\rho(\mathbf{r})d^3r} . \end{equation*}$

What is the meaning of the $\epsilon_i^\mathrm{KS}$ ? Due to the similarities between the Hartree-Fock and Kohn-Sham equations one is inclined to assume that also the Kohn-Sham eigenvalues can be directly related to excitation energies as the ionization potentials or electron affinities. However, there is no full formal justification for this, and one can only show that

for exact DFT, meaning for the exact $E_\mathrm{xc}\mathrm{,}$ or $V_\mathrm{xc}\mathrm{,}$ the highest occupied Kohn-Sham eigenvalue is the negative of the ionization potential
(49) $\begin{equation*} \epsilon_\mathrm{HOMO}^\mathrm{KS} = -\mathrm{IP} \end{equation*}$

and similarly

(50) $\begin{equation*} \epsilon_\mathrm{LUMO}^\mathrm{KS} = -\mathrm{EA} \end{equation*}$
in general, the eigenvalue $\epsilon_i^\mathrm{KS}$ correspond to the change of the total energy with respect to a change of occupation $n_i$ of the $i\mathrm{-th}$ Kohn-Sham orbital (Janak’s theorem)
(51) $\begin{equation*} \epsilon_i^\mathrm{KS} = \frac{\partial E}{\partial N_i} . \end{equation*}$

Comparison of single-particle and excitation energies

Example 1: Carbonmonoxide

Exp. HF DFT-PBE DFT-B3LYP

IP
EA 14.0
1.3 13.1
-3.2 13.8
-2.8 14.2
-2.6

$\epsilon_4$ -19.8 -21.9 -13.9 -15.7

$\epsilon_{5/6}$ -17.0 -17.4 -11.7 -13.1

$\epsilon_\mathrm{HOMO}$ -14.0 -15.0 -8.8 -10.4

$\epsilon_\mathrm{LUMO}$ -1.3 +4.0 -1.7 -0.8

The above table lists results of several Hartree-Fock or Kohn-Sham energies $\epsilon_i$ for a single CO molecule in vacuum compared to experimental energies (PUT REF). One can see that the Hartree-Fock energies for the occupied states agree within 1-2 eV with the experimental reference. So, while Koopman’s theorem allows to make a formal link between the $\epsilon_i$ and the ionization energies, the lack of electron correlation in the effective potential leads to a disagreement with measurements. The agreement of the LUMO and higher energies is usually poor. In the DFT results with both a semi-local (PBE) and hybrid (B3LYP) functional, one finds that both (massively) underestimate the experimental data. For the occupied levels, there is an incomplete cancellation of the repulsive Coulomb self-interaction in the Hartree term by the approximate exchange-correlation potential, yielding too high energies. Also the too repulsive potential can lead to an over-delocalization of the single-particle orbitals.
To obtain more reliable values for the ionization potential and electron affinity from any effective single-particle theory, it is advised not to use the $\epsilon_i$ but determine the total energie differences (this is often referred to as the $\Delta\mathrm{SCF}$ approach):

(52) $\begin{eqnarray*} \mathrm{IP} &=& E[n-1] - E[n]\\ \mathrm{EA} &=& E[n] - E[n+1] \end{eqnarray*}$

The resulting energies are listed at the bottom of the table. Here, we see that the DFT IPs agree very well with the reference value, and the influence of the different functionals seems to be on the order of 0.1eV, roughly an order of magnitude smaller than on the $\epsilon_i$ .
Regarding the electron affinity, all calculations predict a negative electron affinity, i.e., they predict that CO does not want to add an additional electron, which contradicts experimental observations. In fact, the electron affinity of CO is very unique and seems to be a real problem even for higher order quantum-chemical methods.
Example: Fullerene C₆₀

Exp. HF DFT-PBE DFT-B3LYP

IP
EA 7.6
2.7 7.4
1.1 7.5
2.6 7.6
2.5

$\epsilon_\mathrm{HOMO-1}$ -9.0 -9.7 -6.9 -7.6

$\epsilon_\mathrm{HOMO}$ -7.6 -7.9 -5.8 -6.3

$\epsilon_\mathrm{LUMO}$ -2.7 -0.8 -4.1 -3.5

For a considerably bigger molecule, such as C₆₀, we find basically the same trends as in CO with respect to the accuracy of the $\epsilon_i$ . DFT-based $\Delta\mathrm{SCF}$ calculations, on the other hand, agree within 0.2eV with the reference data. The respective Hartree-Fock electron affinity in contrast is off by more than 1eV.

	Exp.	HF	DFT-PBE	DFT-B3LYP
IP EA	14.0 1.3	13.1 -3.2	13.8 -2.8	14.2 -2.6
$\epsilon_4$	-19.8	-21.9	-13.9	-15.7
$\epsilon_{5/6}$	-17.0	-17.4	-11.7	-13.1
$\epsilon_\mathrm{HOMO}$	-14.0	-15.0	-8.8	-10.4
$\epsilon_\mathrm{LUMO}$	-1.3	+4.0	-1.7	-0.8

	Exp.	HF	DFT-PBE	DFT-B3LYP
IP EA	7.6 2.7	7.4 1.1	7.5 2.6	7.6 2.5
$\epsilon_\mathrm{HOMO-1}$	-9.0	-9.7	-6.9	-7.6
$\epsilon_\mathrm{HOMO}$	-7.6	-7.9	-5.8	-6.3
$\epsilon_\mathrm{LUMO}$	-2.7	-0.8	-4.1	-3.5

Lessons:

No method is truly exact!
Hartree-Fock has exact exchange but neglects electron correlations
- energies of occupied orbitals are in reasonable agreement with experimental ionization energies
- unoccupied (virtual) orbital energies are a bad estimate for EAs
DFT is only a formally exact representation of the many-body problem, but we “cheated” by putting all quantum effects into the unknown exchange-correlation functional. Due to the approximate nature of the approximate potentials, we introduce errors, typically:
- - occupied energy levels underestimate the IPs
  - unoccupied energy levels underestimate the EAs
$\Rightarrow$ about $50\%$ underestimate of HOMO-LUMO gaps
HOMO and LUMO energies can be related to electron removal and addition energies. They can be interpreted as representative levels for hole/electron transfer, though a more reliable estimate is obtained by a $\Delta\mathrm{SCF}$ calculation.
The HOMO-LUMO gap (or IP-EA) is not the excitation energy relevant for, e.g., photoabsorption. Such an absorption process is the promotion of one electron from the HOMO to the LUMO (in the most simple case). This is not identical to an independent removal of an electron from the HOMO and addition to the LUMO.Example: C₆₀ from experiment
(53) $\begin{eqnarray*} \mathrm{IP}-\mathrm{EA}&=&\epsilon_\mathrm{LUMO}-\epsilon_\mathrm{HOMO}=4.9\mathrm{eV}\\ \Omega^\mathrm{abs}&=& E^\mathrm{excited}-E_0 = 1.9\mathrm{eV}\\ \Rightarrow \Delta E &=& \mathrm{IP}-\mathrm{EA} - \Omega^\mathrm{abs} = 3.0\mathrm{eV} \end{eqnarray*}$

The difference $\Delta E$ between the free interlevel transition (the HOMO-LUMO gap) and the excitation energy $\Omega^\mathrm{abs}$ is the electron-hole binding energy. An optical excitation, i.e., a coupled two-particle excitation, can never be described by independent effective single-particles.
Apart from over-interpreting the quantitative results of the one has to be careful about
- PES for van-der-Waals bonded systems, e.g., a benzene dimer. Here, HF for instance does not predict any bound dimer configuration (as a function of separation). DFT-PBE does in this case but in general is not reliable;
- molecules with strongly structures electron densities is larger molecular structures, such as donor-acceptor architectures or molecules containing long conjugated $\pi$ systems, due to DFT over-delocalization;
- not introducing new physics/chemistry by using inappropriate computational parameters, like choice of the DFT functional, basis sets, point samplings, convergence criteria etc.

Back to Table of Contents.

I.4 Fermi’s Golden Rule

In the following, we want to consider a two-level system. The two states $\vert 1 \rangle$ and $\vert 2 \rangle$ shall correspond to two eigenstates of the system, in which subsystems 1 and 2 do not interact (diabatic states). Assuming that the system is initially prepared in state $\vert 1\rangle$ , we would like to know the rate of transfer into state $\vert 2 \rangle$ if we now switch on a (small) interaction.
First, we can express the Hamiltonian of the full system in the basis of the two diabatic states

(54) $\begin{eqnarray*} \hat{H} &=& \langle 1 \vert \hat{H} \vert 1 \rangle \vert 1 \rangle\langle 1 \vert + \langle 2 \vert \hat{H} \vert 2 \rangle \vert 2 \rangle\langle 2 \vert + \langle 1 \vert \hat{H} \vert 2 \rangle \vert 1 \rangle\langle 2 \vert + \langle 2 \vert \hat{H} \vert 1 \rangle \vert 2 \rangle\langle 1 \vert \\ &=& E_1 \rangle \vert 1 \rangle\langle 1 \vert + E_2 \rangle \vert 2 \rangle\langle 2 \vert + V_{12} \rangle \vert 1 \rangle\langle 2 \vert + V_{21} \rangle \vert 2 \rangle\langle 1 \vert , \end{eqnarray*}$

where $E_1$ and $E_2$ are the energies of the two independent subsystems and $V_{12}=V_{21}^*$ is the coupling element between them. The fact that the system is supposed to be initially in state $\vert 1 \rangle$ can be described in terms of the occupation probability $P_1(t=0)=1$ . Evidently, $P_2(t=0)=0$ . The question now is how this occupation probability changes with time? We can associate the change of occupation probability with a transfer rate $k_{1\rightarrow 2}$ :

(55) $\begin{equation*} \frac{\partial}{\partial t} P_1(t) = -k_{1\rightarrow 2} P_1(t) . \end{equation*}$

The time evolution of the initial state can be described using the propagator

(56) $\begin{equation*} \hat{U}(t)=\exp{\left(-\frac{i}{\hbar}\hat{H}t\right)} \end{equation*}$

according to

(57) $\begin{equation*} \vert \Psi(t)\rangle = \hat{U}(t) \vert 1 \rangle . \end{equation*}$

To proceed, we form matrix elements of $\hat{U}(t)$ with the diabatic states $\vert \nu \rangle$ and $\vert \mu \rangle$

(58) $\begin{equation*} A_{\nu\mu}(t) = \theta(t) \langle \nu \vert \hat{U}(t) \vert \mu \rangle . \end{equation*}$

These transition amplitudes are a measure of how much of the initial state $\vert \nu \rangle$ is still contained in the propagated state $\hat{U}(t)\vert\mu\rangle$ . The time derivative of the transition amplitudes read

(59) $\begin{eqnarray*} \frac{\partial}{\partial t} A_{\nu\mu} &=& \delta(t)\delta_{\nu\mu}\delta_{\mu 1} - \frac{i}{\hbar}\theta(t) \langle \nu \vert \hat{H}\hat{U}(t) \vert \mu \rangle\\ &=& \delta(t)\delta_{\nu\mu}\delta_{\mu 1} - \frac{i}{\hbar}\theta(t) \sum_k{\langle \nu \vert \hat{H} \vert k\rangle\langle k \vert \hat{U}(t)\vert\mu\rangle} \end{eqnarray*}$

which leads to

(60) $\begin{equation*} i\hbar \frac{\partial}{\partial t} A_{\nu\mu} = i \hbar \delta(t)\delta_{\nu\mu}\delta_{\mu 1} + \sum_k{\langle \nu \vert \hat{H} \vert k\rangle A_{k\mu} } . \end{equation*}$

Fourier transforming the last equation to frequency space, we obtain

(61) $\begin{equation*} \hbar \omega A_{\mu\nu}(\omega) = i\hbar \delta_{\nu\mu} \delta_{\mu 1} + \sum_k{\langle \nu \vert \hat{H} \vert k\rangle A_{k\mu}(\omega)} . \end{equation*}$

Now, let us consider the special case of $\nu=\mu=1\mathrm{,}$ i.e., the survival amplitude of state $\vert 1 \rangle\mathrm{:}$

(62) $\begin{equation*} \hbar \omega A_{11}(\omega) = i\hbar + E_1 A_{11}(\omega) + V_{12}A_{21}(\omega) \end{equation*}$

for which we obviously need to know $A_{21}(\omega)\mathrm{:}$

(63) $\begin{eqnarray*} \hbar \omega A_{21}(\omega) &=& V_{21}A_{11}(\omega) + E_2 A_{21}(\omega)\\ \Rightarrow A_{21}(\omega) &=& \frac{V_{21}}{\hbar \omega -E_2 + i\epsilon} \end{eqnarray*}$

where we introduced an infinitesimal imaginary part $i\epsilon$ to allow inversion of the equation. We will take the limit $\epsilon \to 0$ later. So, finally, we have found a closed expression for the survival amplitude

(64) $\begin{equation*} A_{11}(\omega) = i\hbar \left\{ \hbar \omega - E_1 + \frac{\vert V_{12} \vert^2}{\hbar\omega - E_2 + i \epsilon}\right\}^{-1} \\ \end{equation*}$

Now, consider

(65) $\begin{eqnarray*} \lim_{\epsilon \to 0}{\left[\frac{\vert V_{12} \vert^2}{\hbar\omega - E_2 + i \epsilon} \right]} &=& \lim_{\epsilon \to 0}{\left[ \frac{\vert V_{12}\vert^2 (\hbar \omega - E_2)}{(\hbar \omega -E_2)+\epsilon^2} - i \frac{\epsilon \vert V_{12}\vert^2}{(\hbar\omega - E_2)+\epsilon^2} \right]} \\ &=& P\left( \frac{\vert V_{12}\vert^2}{\hbar \omega -E_2}\right) - i\pi \vert V_{12}\vert^2 \delta(\hbar \omega - E_2)\\ &=& \hbar \Omega(\omega) - i\hbar \Gamma(\omega), \end{eqnarray*}$

where $P(\ldots)$ denotes the Principal value. The time-dependent occupation probability of the initial state is then determined as

(66) $\begin{equation*} P_1(t) = \vert A_{11}(t)\vert^2 = \left\vert \int{e^{-i\omega t}\frac{i\hbar}{\hbar \omega - E_1 + \hbar \Omega(\omega)-i\hbar\Gamma(\omega)}\frac{d\omega}{2\pi}} \right\vert^2 . \end{equation*}$

To evaluate this integral, we make the following assumptions:

The variation of $\Omega(\omega)$ and $\Gamma(\omega)$ is weak for $\omega \approx E_1/\hbar$ .
The $\omega$ -dependence results mainly from the resonance in $\omega=E_1/\hbar$ .

Then, we can use the residue theorem to finally obtain

(67) $\begin{equation*} P_1(t) \approx \left\vert \int{e^{-i\omega t}\frac{i\hbar}{\hbar \omega - E_1 + \hbar \Omega\left(E_1/\hbar\right)-i\hbar\Gamma\left(E_1/\hbar\right)}\frac{d\omega}{2\pi}} \right\vert^2 = \theta(t)\exp{\left[-2\Gamma\left(E_1/\hbar\right)t\right]} \end{equation*}$

Taking the time derivative

(68) $\begin{equation*} \frac{\partial}{\partial t} P_1(t) = -2\Gamma\left(E_1/\hbar\right) P_1(t) \end{equation*}$

and comparing this to Eq.55 one can read off the transfer rate

(69) $\begin{equation*} k_{1\to 2} = \frac{2\pi}{\hbar} \vert V_{12}\vert^2 \delta(E_1-E_2) , \end{equation*}$

in which the term $\delta(E_1-E_2)$ ensures energy conservation. For a set of final states $\vert \alpha\rangle$ instead of a single final states, the rate generalizes to Fermi’s Golden Rule

(70) $\begin{equation*} k = \frac{2\pi}{\hbar} \sum_\alpha{\vert V_{1\alpha}\vert^2 \delta(E_1-E_\alpha)} . \end{equation*}$

Back to Table of Contents.

I.5 Harmonic approximation for nuclear motion

Assume that we have solved the electronic problem for different nuclear conformations, i.e., we have determined the potential energy surface $U_\nu(\mathbf{R})$ (or the effective potential) for the nuclear dynamics. If the molecule exists, there is a stable equilibrium configuration $\mathbf{R}_0$ with $V_0=U_0(\mathbf{R}_0)$ and we can express the dynamical nuclear coordinates using deviations $\mathbf{u}$ from this configuration

(71) $\begin{equation*} \mathbf{R} = \mathbf{R}_0 + \mathbf{u} . \end{equation*}$

We can then Taylor expand the effective potential around $\mathbf{R}_0$ (nucleus $\alpha$ , Cartesian coordinate $l$ )

(72) $\begin{equation*} U(\mathbf{R}) = \underbrace{U(\mathbf{R}_0)}_{=V_0} + \underbrace{ \sum_{\alpha l}{\left. \frac{\partial U(\mathbf{R})}{\partial R_{\alpha l}}\right\vert_{\mathbf{R}_0}}u_{\alpha l}}_{=0 \mathrm{,\,stationary point}}} + \frac{1}{2} \sum_{\alpha l, \beta m}{\underbrace{\left.\frac{\partial^2 U(\mathbf{R})}{\partial R_{\alpha l}\partial R_{\beta m}}\right\vert_{\mathbf{R}_0}}_{=\kappa_{\alpha l,\beta m}\mathrm{, Hessian}} u_{\alpha l} u_{\beta m}} + \ldots \end{equation*}$

In this harmonic approximation, the effective potential for the nuclear dynamics is given by

(73) $\begin{equation*} U(\mathbf{R}) = V_0 + \frac{1}{2} \sum_{\alpha l,\beta m}{\kappa_{\alpha l,\beta m} u_{\alpha l} u_{\beta m}} \end{equation*}$

and consequently the nuclear Hamiltonian by

(74) $\begin{equation*} \hat{H}_\mathrm{nuc} = \sum_{\alpha l}{\frac{\hat{P}_{\alpha l}^2}{2M_\alpha}} + V_0 + \frac{1}{2} \sum_{\alpha l,\beta m}{\kappa_{\alpha l,\beta m} u_{\alpha l} u_{\beta m}} . \end{equation*}$

Introducing mass-weighted normal coordinates

(75) $\begin{eqnarray*} \bar{u}_{\alpha l} &=& \sqrt{M_\alpha} u_{\alpha l} \\ \bar{P}_{\alpha l} &=& \frac{1}{\sqrt{M_\alpha}}P_{\alpha l} \end{eqnarray**}$

one can rewrite the Hamiltonian as

(76) $\begin{equation*} \hat{H}_\mathrm{nuc} = \frac{1}{2} \left(\mathbf{\bar{P}}^\dagger\mathbf{\bar{P}} + \mathbf{\bar{u}}^\dagger \mathbf{\underline{D}} \mathbf{\bar{u}} \right) \end{equation*}$

where

(77) $\begin{equation*} D_{\alpha l, \beta m} = \frac{1}{\sqrt{M_\alpha M_\beta}} \kappa_{\alpha l, \beta m} \end{equation*}$

is the dynamical matrix. This matrix can be diagonalized by an orthogonal transformation

(78) $\begin{equation*} \mathbf{\underline{\Omega}} = \mathbf{\underline{C}} \mathbf{\underline{D}} \mathbf{\underline{C}}^\dagger = \left(\begin{array}{ccc} \omega_1^2 & 0 & 0 \\ 0 & \ddots & 0 \\ 0 & 0 & \omega_{3N}^2 \end{array} \right) \end{equation*}$

in which the $3N$ $\omega_i$ are normal-mode frequencies. Applying the transformation $\mathbf{\underline{C}}$ to the coordinates and impulses

(79) $\begin{eqnarray*} \mathbf{\tilde{u}} &=& \mathbf{\underline{C}} \mathbf{\bar{u}}\\ \mathbf{\tilde{P}} &=& \mathbf{\underline{C}} \mathbf{\bar{P}} \end{eqnarray*}$

finally yields the nuclear Hamiltonian in the form

(80) $\begin{equation*} \hat{H}_\mathrm{nuc} = \frac{1}{2} \sum_{i=1}^{3N}{\left( \tilde{P}_i^2 + \omega_i^2 \tilde{u}_i^2\right)} . \end{equation*}$

This Hamiltonian is evidently a sum over $3N$ independent harmonic oscillators, for each of which the solution is known, i.e., $\chi_{\lambda_i}(\tilde{u}_i)$ is the $\lambda_i\mathrm{-th}$ eigenfunction of the harmonic oscillator and the total nuclear wave function is given by