“The theory of fermions is one of the great triumphs of twentieth-century physics. Most of the credit belongs to Paul Dirac, who started, like Einstein, with some simple assumptions and laid the foundations for the Pauli exclusion principle of chemistry, Fermi statistics in solids, and antimatter in particle physics.”

Deduce the Dirac equation

The basic problem of Schrodinger equation is it doesn’t contain relativistic effect, since if we wrote: \(\begin{align*} i\hbar \frac{\partial}{\partial t}\psi=H\psi=\left[\frac{-\hbar^2\nabla^2}{2m}+U\right]\psi \end{align*}\) we can see the kinetic terms is not relativistically invariant. The conserved quantities in special relativity is $E^2=(mc^2)^2+(cp)^2$. However, if the time derivitive term has second order. Then the equation will become like: \(\begin{align} \left(i\hbar\frac{\partial}{\partial t}\right)^2\psi&=-\hbar^2\frac{\partial^2}{\partial t^2}\psi\\ &=\left[m^2c^4+c^2(-\hbar^2\nabla^2)\right]\psi\\ &=E^2\psi\\ \Leftrightarrow -\frac{1}{c^2}\frac{\partial^2}{\partial t^2}\psi&=\frac{m^2c^2}{\hbar^2}\psi-\nabla^2\psi \end{align}\) There will be two problem appears:

1. Since this equation is second order, the zero and first order term of $\psi$ can be designed arbitrarily, therefore, the $\psi$ lost it’s physical intepretation where it’s amplitude equals the probability.

2. The solution of this equation will looks like $\psi(x,t)=\psi_0\cos(\omega t)\cos(kx)$, which means particle can disappear if $\cos(\omega t)=0$.

So “Dirac believed strongly the equation should be dependent in time linearly.” (btw the above quadratic equation is known as Klein-Gordon equation.)

Dirac divided \(E^2-(mc^2)^2-|c\vec{p}|^2=(E+\alpha_0mc^2+c\vec{\alpha}\cdot\vec{p})(E-\alpha_0mc^2-c\vec{\alpha}\cdot\vec{p})=0\) which can be done if ${\alpha_i,\alpha_j}=2\delta_{ij}$ is anticommutational. So Paili matrices is a natrual choice. Then we have Dirac equation: \(i\hbar\frac{\partial}{\partial t}|\psi\rangle=H|\psi\rangle=(\alpha_0mc^2+c\vec{\alpha}\cdot\vec{p})|\psi\rangle\) where \(\alpha_0= \left( \begin{matrix} E & 0 \\ 0 & -E \end{matrix} \right), \quad \alpha_i= \left( \begin{matrix} 0 & \sigma_i \\ \sigma_i & 0 \end{matrix} \right)\) That $\sigma_i$ is Pauli spin matrices. This equation will have both positive an negative energy solution, which indicates positrons in the Dirac sea.

The Dirac equation and spin energies

Before we only includes the kinetic term in Hamiltonian, now we can add in the magnetic term and the electrostatic term, to see what happen. Aware the conservative energy \(E^2=(mc^2)^2-|c\vec{\alpha}\cdot\vec{p}|^2.\)

Magnetic Effect To account for magnetic field, we can have a vector potential, and then: $E^2=(mc^2)^2+|c\vec{\alpha}\cdot(\vec{p}-q\vec{A})|^2$. Assume $|c\vec{\alpha}\cdot(\vec{p}-q\vec{A})|$ is comparably smaller than the first, we have $E=\sqrt{(mc^2)^2+(cp)^2}\approx mc^2+p^2/2m$, by $\alpha_x\alpha_y=i\sigma_z$ we have: \(\begin{align*} E^2&=(mc^2)^2+c^2|\vec{\alpha}\cdot(\vec{p}-q\vec{A})|^2\\&=(mc^2)^2+c^2\left|\sum_{i}\alpha_i(p_i-qA_i)\right|^2\\&=(mc^2)^2+c^2\left[\sum_{i}\alpha_i^2(p_i-qA_i)^2+\sum_{i\neq j}\alpha_i\alpha_j(p_i-qA_i)(p_j-qA_j)\right]\\&=(mc^2)^2+c^2|\vec{p}-q\vec{A}|^2+\hbar c^2\vec{\sigma}\cdot(\nabla\times\vec{A}) \end{align*}\)

Using the first order approximation, We have then:

The first two terms correspond to the standard nonrelativistic Hamiltonian for a particle in a magnetic field, while the last term is a new term that arises from the fact that we needed to introduce the α matrices for the relativistic wave equation. This term also $\frac{q\hbar}{2m}\vec{\sigma}\cdot\vec{B}$ been denoted as the Zeeman effect. The $\frac{q\hbar}{2m}$ are acting as the magnetic momentum. Since magnetic potential is often written as $U_B=\vec{m}\cdot \vec{B}$.

Electrostatic Effect Then we include a electrostatic field instead of the $\vec{A}$. We have \(E|\psi\rangle=\left(\alpha_0mc^2+c\vec{\alpha}\cdot\vec{p}+U(\vec{r})\right)|\psi\rangle\)

By definition of $\alpha_0$ and $\alpha$, the equation above can be split into two, which is: \(\begin{align} E|\psi_1\rangle&=mc^2|\psi_1\rangle+c\vec{\sigma}\cdot\vec{p}|\psi_2\rangle+U(\vec{r})|\psi_1\rangle\\ E|\psi_2\rangle&=-mc^2|\psi_2\rangle+c\vec{\sigma}\cdot\vec{p}|\psi_1\rangle+U(\vec{r})|\psi_2\rangle \end{align}\) Let $E’=E-mc^2$, so: \(\begin{align} E'|\psi_1\rangle&=c\vec{\sigma}\cdot\vec{p}|\psi_2\rangle+U(\vec{r})|\psi_1\rangle\\ E'|\psi_2\rangle&=c\vec{\sigma}\cdot\vec{p}|\psi_1\rangle+\left[U(\vec{r})-2mc^2\right]|\psi_2\rangle \end{align}\)

Noticing that all the relativistically effect of $\psi_1$ is included via coupling of $\psi_2$.

Rewrite (8) and approximate with second term in tayler series, we have: \(|\psi_2\rangle=\frac{1}{2mc^2+E'-U}c\vec{\sigma}\cdot\vec{p}|\psi_1\rangle\approx c\left[\frac{1}{2mc^2}-\frac{E'-U}{(2mc^2)^2}\right]\vec{\sigma}\cdot\vec{p}|\psi_1\rangle=\frac{1}{2mc}\left(1-\frac{E'}{2mc^2}+\frac{U}{2mc^2}\right)\vec{\sigma}\cdot\vec{p}|\psi_1\rangle\)

substitude in (7) to obtain: \(\begin{align*} E'&=\frac{1}{2mc}\vec{\sigma}\cdot\vec{p}\left(1-\frac{E'}{2mc^2}+\frac{U}{2mc^2}\right)\vec{\sigma}\cdot\vec{p}+U(\vec{r})\\ &=\frac{1}{2mc}(\vec{\sigma}\cdot\vec{p})(\vec{\sigma}\cdot\vec{p})-\frac{E'}{4m^2c^4}(\vec{\sigma}\cdot\vec{p})(\vec{\sigma}\cdot\vec{p})+(\vec{\sigma}\cdot\vec{p})\frac{U}{4m^2c^4}(\vec{\sigma}\cdot\vec{p})+U(\vec{r})\\ &=\frac{p^2}{2m}\left(1-\frac{E'-U}{2mc^2}\right)+U(\vec{r})-\frac{1}{4m^2c^2}(i\hbar\nabla U)\cdot\vec{p}+\frac{\hbar}{4m^2c^2}\vec{\sigma}\cdot(\nabla U\times\vec{p}). \end{align*}\)

The first term gives the kinetic energy with a relativistic correction for the mass; the third term is relativistic correction for the potential energy. The last term is the standard spinorbit energy.

Spin-orbital Coupling

When considering the central potential $U(r)$, where $\nabla U=\partial_rU\hat{r}$, the last term can be simplified as:

Where $\vec{L}=\vec{r}\times\vec{p}$, $\vec{S}=\frac{\hbar}{2}\vec{\sigma}$. This formula appears very often in practical calculation in DFT, since the zeros order approximation to the SOC term also gives similar format.

Spin-orbital Coupling in DFT

In various DFT implementations, SOC are treated differently. The most direct way is to solve the four component Dirac equation. This requires no approximation, but need to include four times the size of the basis. This is implemented in the method mostly based on all electron DFT. Another way is use the approximated equations that have two fold wave function as described above.(Also in all-electron calculations mostly):

We can see that, the correction to the kinetical energy and potential does not contain spin, therefore, including only those corrections can be easily incooperated in any code, without considering the spin nondegeneracy.

Further simplification is possible utilizing the pesudopotentials (PP). Since the SOC term is dominated by the gradient of the potential. Therefore, the core electrons have larger split comparing to the valance ones. By discarding the valance’s SOC contribution, we can include all the SOC contribution in the PPs. The decomposition of effect can be justfies with:

\(\frac{\hbar}{4m^2c^2}\vec{\sigma}\cdot(\nabla U\times\vec{p})\approx \sum_i\xi_i(r)\vec{L}_i\cdot\vec{S}_i\) Where $\xi_i(r)=\frac{1}{2m^2c^2}\frac{\partial U}{\partial r}$ is often small outside the atomic sphere, while inside, the potential $U(r)$ is approximately spherical, therefore approximate as spherical.

The PP generated from the fully dirac equation can be decomposed into the scalar part and the angular dependent part, where the SOC effect are contains in the angular dependent part. Using only the scalar correction in of relativistic effect of the kinetical energy and potentials are denoted as scalar-relativistic (SR). The further inlusion of the SOC term (and possibly other method to generate the PP, like solving dirac equations directly, or 2-component formalisms), is denoted as the fully-relativistic (FR).

Details of how to incooperate FR PP and the SOC Hamiltonian in DFT can be find in [2] briefly, the pseudopotential is decomposed as:

\(\begin{align*} V^{FR}&=V^{SR}+V^{SO}\\ &=\sum_{lm}|v_{lm}^{SR}\rangle\langle v_{lm}^{SR}|+\sum_{lm}\left[\frac{1}{4}l(l+1)-\frac{1}{2}L\cdot S\right]|v_{lm}^{SO}\rangle\langle v_{lm}^{SO}| \end{align*}\) Where $v_{SO}$ and $v_{SR}$ are transformed from the nonlocal Kleinman–Bylander [3] form of the fully relativistic pseudopotential of $v_{lJ}$.

Further more, the calculation of the SOC can even further simplified by discard the off diagonal SOC matrix elements that belong to different atoms. The result is named as the onsite approximation as in [4].

Reference ###:

[1] Lenthe, E. van, Evert-Jan Baerends, and Jaap G. Snijders. “Relativistic regular two‐component Hamiltonians.” The Journal of chemical physics 99.6 (1993): 4597-4610.

[2] Cuadrado, R., and J. I. Cerdá. “Fully relativistic pseudopotential formalism under an atomic orbital basis: spin–orbit splittings and magnetic anisotropies.” Journal of Physics: Condensed Matter 24.8 (2012): 086005.

[3] Kleinman L and Bylander D M 1982 Phys. Rev. Lett. 48 1425

[4] Cuadrado, Ramón, et al. “Validity of the on-site spin-orbit coupling approximation.” Physical Review B 104.19 (2021): 195104.