X-Git-Url: https://hackdaworld.org/gitweb/?p=lectures%2Flatex.git;a=blobdiff_plain;f=physics_compact%2Fsolid.tex;h=5e69a799a18a5174921aabe05797cc83202675e6;hp=20ecfef44e2132a8f2f8eb73c9dd44480b5cf773;hb=6a26b9f5593acc0bf19241b2fe79f1acf51fb03e;hpb=a9f93985e52272cdecf902c4b559173a80f4b41d

diff --git a/physics_compact/solid.tex b/physics_compact/solid.tex
index 20ecfef..5e69a79 100644
--- a/physics_compact/solid.tex
+++ b/physics_compact/solid.tex
@@ -2,6 +2,10 @@
 
 \chapter{Atomic structure}
 
+\chapter{Reciprocal lattice}
+
+Example of primitive cell ...
+
 \chapter{Electronic structure}
 
 \section{Noninteracting electrons}
@@ -18,30 +22,441 @@
 
 \subsubsection{Hohenberg-Kohn theorem}
 
-Considering a system with a nondegenerate ground state, there is obviously only one ground-state charge density $n_0(\vec{r})$ that correpsonds to a given potential $V(\vec{r})$.
+The Hamiltonian of a many-electron problem has the form
+\begin{equation}
+H=T+V+U\text{ ,}
+\end{equation}
+where
+\begin{eqnarray}
+T & = & \langle\Psi|\sum_{i=1}^N\frac{-\hbar^2}{2m}\nabla_i^2|\Psi\rangle\\
+  & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} d\vec{r}' \,
+        \langle \Psi | \vec{r} \rangle \langle \vec{r} |
+        \nabla_i^2
+        | \vec{r}' \rangle \langle \vec{r}' | \Psi \rangle\\
+  & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} d\vec{r}' \,
+        \langle \Psi | \vec{r} \rangle \nabla_{\vec{r}_i}
+        \langle \vec{r} | \vec{r}' \rangle
+        \nabla_{\vec{r}'_i} \langle \vec{r}' | \Psi \rangle\\
+  & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} d\vec{r}' \,
+        \nabla_{\vec{r}_i} \langle \Psi | \vec{r} \rangle
+        \delta_{\vec{r}\vec{r}'}
+        \nabla_{\vec{r}'_i} \langle \vec{r}' | \Psi \rangle\\
+  & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} \,
+        \nabla_{\vec{r}_i} \Psi^*(\vec{r}) \nabla_{\vec{r}_i} \Psi(\vec{r})
+        \text{ ,} \\
+V & = & \int V(\vec{r})\Psi^*(\vec{r})\Psi(\vec{r})d\vec{r} \text{ ,} \\
+U & = & \frac{1}{2}\int\frac{1}{\left|\vec{r}-\vec{r}'\right|}
+        \Psi^*(\vec{r})\Psi^*(\vec{r}')\Psi(\vec{r}')\Psi(\vec{r})
+        d\vec{r}d\vec{r}'
+\end{eqnarray}
+represent the kinetic energy, the energy due to the external potential and the energy due to the mutual Coulomb repulsion.
+
+\begin{remark}
+As can be seen from the above, two many-electron systems can only differ in the external potential and the number of electrons.
+The number of electrons is determined by the electron density.
+\begin{equation}
+N=\int n(\vec{r})d\vec{r}
+\end{equation}
+Now, if the external potential is additionally determined by the electron density, the density completely determines the many-body problem.
+\end{remark}
+
+Considering a system with a nondegenerate ground state, there is obviously only one ground-state charge density $n_0(\vec{r})$ that corresponds to a given potential $V(\vec{r})$.
+\begin{equation}
+n_0(\vec{r})=\int \Psi_0^*(\vec{r},\vec{r}_2,\vec{r}_3,\ldots,\vec{r}_N)
+                  \Psi_0(\vec{r},\vec{r}_2,\vec{r}_3,\ldots,\vec{r}_N)
+             d\vec{r}_2d\vec{r}_3\ldots d\vec{r}_N
+\end{equation}
 In 1964, Hohenberg and Kohn showed the opposite and far less obvious result \cite{hohenberg64}.
+
+\begin{theorem}[Hohenberg / Kohn]
 For a nondegenerate ground state, the ground-state charge density uniquely determines the external potential in which the electrons reside.
-The proof presented by Hohenberg and Kohn proceeds by {\em reductio ad absurdum}.
+\end{theorem}
 
+\begin{proof}
+The proof presented by Hohenberg and Kohn proceeds by {\em reductio ad absurdum}.
 Suppose two potentials $V_1$ and $V_2$ exist, which yield the same electron density $n(\vec{r})$.
 The corresponding Hamiltonians are denoted $H_1$ and $H_2$ with the respective ground-state wavefunctions $\Psi_1$ and $\Psi_2$ and eigenvalues $E_1$ and $E_2$.
 Then, due to the variational principle (see \ref{sec:var_meth}), one can write
 \begin{equation}
-E_1=\langle \Psi_1 | H_1 | \Psi_1 \rangle < \langle \Psi_2 | H_1 | \Psi_2 \rangle
+E_1=\langle \Psi_1 | H_1 | \Psi_1 \rangle <
+\langle \Psi_2 | H_1 | \Psi_2 \rangle \text{ .}
+\label{subsub:hk01}
 \end{equation}
-Expressing $H_1$ by $H_2+H_1-H_2$
+Expressing $H_1$ by $H_2+H_1-H_2$, the last part of \eqref{subsub:hk01} can be rewritten:
 \begin{equation}
 \langle \Psi_2 | H_1 | \Psi_2 \rangle =
 \langle \Psi_2 | H_2 | \Psi_2 \rangle +
 \langle \Psi_2 | H_1 -H_2 | \Psi_2 \rangle
 \end{equation}
-and the fact that the two Hamiltonians, which describe the same number of electrons, differ only in the potential
+The two Hamiltonians, which describe the same number of electrons, differ only in the potential
 \begin{equation}
 H_1-H_2=V_1(\vec{r})-V_2(\vec{r})
 \end{equation}
-one obtains
+and, thus
 \begin{equation}
 E_1<E2+\int n(\vec{r}) \left( V_1(\vec{r})-V_2(\vec{r}) \right) d\vec{r}
+\text{ .}
+\label{subsub:hk02}
+\end{equation}
+By switching the indices of \eqref{subsub:hk02} and adding the resulting equation to \eqref{subsub:hk02}, the contradiction
+\begin{equation}
+E_1 + E_2 < E_2 + E_1 +
+\underbrace{
+\int n(\vec{r}) \left( V_1(\vec{r})-V_2(\vec{r}) \right) d\vec{r} +
+\int n(\vec{r}) \left( V_2(\vec{r})-V_1(\vec{r}) \right) d\vec{r}
+}_{=0}
+\end{equation}
+is revealed, which proofs the Hohenberg Kohn theorem.% \qed
+\end{proof}
+
+\section{Electron-ion interaction}
+
+\subsection{Pseudopotential theory}
+
+The basic idea of pseudopotential theory is to only describe the valence electrons, which are responsible for the chemical bonding as well as the electronic properties for the most part.
+
+\subsubsection{Orthogonalized planewave method}
+
+Due to the orthogonality of the core and valence wavefunctions, the latter exhibit strong oscillations within the core region of the atom.
+This requires a large amount of planewaves $\ket{k}$ to adequatley model the valence electrons.
+
+In a very general approach, the orthogonalized planewave (OPW) method introduces a new basis set
+\begin{equation}
+\ket{k}_{\text{OPW}} = \ket{k} - \sum_t \ket{t}\bra{t}k\rangle \text{ ,}
+\end{equation} 
+with $\ket{t}$ being the eigenstates of the core electrons.
+The new basis is orthogonal to the core states $\ket{t}$.
+\begin{equation}
+\braket{t}{k}_{\text{OPW}} =
+\braket{t}{k} - \sum_{t'} \braket{t}{t'}\braket{t'}{k} =
+\braket{t}{k} - \braket{t}{k}=0
+\end{equation}
+The number of planewaves required for reasonably converged electronic structure calculations is tremendously reduced by utilizing the OPW basis set.
+
+\subsubsection{Pseudopotential method}
+
+Following the idea of orthogonalized planewaves leads to the pseudopotential idea, which --- in describing only the valence electrons --- effectively removes an undesriable subspace from the investigated problem.
+
+Let $\ket{\Psi_\text{V}}$ be the wavefunction of a valence electron with the Schr\"odinger equation 
+\begin{equation}
+H \ket{\Psi_\text{V}} = \left(\frac{1}{2m}p^2+V\right)\ket{\Psi_\text{V}}=
+E\ket{\Psi_\text{V}} \text{ .}
+\end{equation}
+\ldots projection operatore $P_\text{C}$ \ldots
+
+\subsubsection{Semilocal form of the pseudopotential}
+
+Ionic potentials, which are spherically symmteric, suggest to treat each angular momentum $l,m$ separately leading to $l$-dependent non-local (NL) model potentials $V_l(r)$ and a total potential
+\begin{equation}
+V=\sum_{l,m}\ket{lm}V_l(\vec{r})\bra{lm} \text{ .}
+\end{equation}
+In fact, applied to a function, the potential turns out to be non-local in the angular coordinates but local in the radial variable, which suggests to call it asemilocal (SL) potential.
+
+Problem of semilocal potantials become valid once matrix elements need to be computed.
+Integral with respect to the radial component needs to be evaluated for each planewave combination, i.e.\  $N(N-1)/2$ integrals.
+\begin{equation}
+\bra{k+G}V\ket{k+G'} = \ldots
+\end{equation}
+
+A local potential can always be separated from the potential \ldots
+\begin{equation}
+V=\ldots=V_{\text{local}}(\vec{r})+\ldots
+\end{equation}
+
+\subsubsection{Norm conserving pseudopotentials}
+
+HSC potential \ldots
+
+\subsubsection{Fully separable form of the pseudopotential}
+
+KB transformation \ldots
+
+\subsection{Spin-orbit interaction}
+
+Relativistic effects can be incorporated in the normconserving pseudopotential method up to but not including terms of order $\alpha^2$ \cite{kleinman80,bachelet82} with $\alpha$ being the fine structure constant.
+This is advantageous since \ldots
+With the solutions of the all-electron Dirac equations, the new pseudopotential reads
+\begin{equation}
+V(\vec{r})=\sum_{l,m}\left[
+\ket{l+\frac{1}{2},m+{\frac{1}{2}}}V_{l,l+\frac{1}{2}}(\vec{r})
+\bra{l+\frac{1}{2},m+{\frac{1}{2}}} +
+\ket{l-\frac{1}{2},m-{\frac{1}{2}}}V_{l,l-\frac{1}{2}}(\vec{r})
+\bra{l-\frac{1}{2},m-{\frac{1}{2}}}
+\right] \text{ .}
+\label{eq:solid:so_bs1}
+\end{equation}
+By defining an averaged potential weighted by the different $j$ degeneracies of the $\ket{l\pm\frac{1}{2}}$ states
+\begin{equation}
+\bar{V}_l(r)=\frac{1}{2l+1}\left(
+l V_{l,l-\frac{1}{2}}(\vec{r})+(l+1)V_{l,l+\frac{1}{2}}(\vec{r})\right)
+\end{equation}
+and a potential describing the difference in the potential with respect to the spin
+\begin{equation}
+V^{\text{SO}}_l(\vec{r})=\frac{2}{2l+1}\left(
+V_{l,l+\frac{1}{2}}(\vec{r})-V_{l,l-\frac{1}{2}}(\vec{r})\right)
+\end{equation}
+the total potential can be expressed as
+\begin{equation}
+V(\vec{r})=\sum_l
+\ket{l,m}\left[\bar{V}_l(\vec{r})+V^{\text{SO}}_l(\vec{r})LS\right]\bra{l,m}
+\text{ ,}
+\label{eq:solid:so_bs2}
+\end{equation}
+where the first term correpsonds to the mass velocity and Darwin relativistic corrections and the latter is associated with the spin-orbit (SO) coupling.
+\begin{proof}
+This can be shown by rewriting the $LS$ operator
+\begin{equation}
+J=L+S \Leftrightarrow J^2=L^2+S^2+2LS \Leftrightarrow
+LS=\frac{1}{2}\left(J^2-L^2-S^2\right)
+\end{equation}
+and corresponding eigenvalue
+\begin{eqnarray}
+j(j+1)-l(l+1)-s(s+1)&=&
+(l\pm\frac{1}{2})(l\pm\frac{1}{2}+1)-l^2-l-\frac{3}{4} \nonumber\\
+&=&
+l^2\pm\frac{l}{2}+l\pm\frac{l}{2}+\frac{1}{4}\pm\frac{1}{2}-l^2-l-\frac{3}{4}
+\nonumber\\
+&=&\pm(l+\frac{1}{2})-\frac{1}{2}=\left\{\begin{array}{rl}
+l & \text{for } j=l+\frac{1}{2}\\
+-(l+1) & \text{for } j=l-\frac{1}{2}
+\end{array}\right.
+\text{ ,}
+\end{eqnarray}
+which, if used in equation~\eqref{eq:solid:so_bs2}, gives the same (diagonal) matrix elements
+\begin{eqnarray}
+\bra{l\pm\frac{1}{2},m\pm\frac{1}{2}}V(\vec{r})
+\ket{l\pm\frac{1}{2},m\pm\frac{1}{2}}&=&
+\bar{V}_l(\vec{r})+V^{\text{SO}}_l(\vec{r})
+\frac{1}{2}\left(l(l+1)-j(j+1)-\frac{3}{4}\right) \nonumber\\
+&=&\bar{V}_l(\vec{r})+\frac{1}{2}V^{\text{SO}}_l(\vec{r})
+\cdot\left\{\begin{array}{cl}
+l & \text{for } j=l+\frac{1}{2}\\
+-(l+1) & \text{for } j=l-\frac{1}{2}
+\end{array}\right. \nonumber\\
+&=&\frac{1}{2l+1}\left(lV_{l,l-\frac{1}{2}}(\vec{r})+
+                       (l+1)V_{l,l+\frac{1}{2}}(\vec{r})\right)+\nonumber\\
+&&\frac{1}{2l+1}
+\left(V_{l,l+\frac{1}{2}}(\vec{r})-V_{l,l-\frac{1}{2}}(\vec{r})\right)
+\cdot\left\{\begin{array}{c}
+l \\
+-(l+1)
+\end{array}\right. \nonumber\\
+&=&\left\{\begin{array}{cl}
+V_{l,l+\frac{1}{2}}(\vec{r}) & \text{for } j=l+\frac{1}{2}\\
+V_{l,l-\frac{1}{2}}(\vec{r}) & \text{for } j=l-\frac{1}{2}
+\end{array}\right.
+\end{eqnarray} 
+as expected and --- in fact --- obtained from equation~\eqref{eq:solid:so_bs1}.
+\end{proof}
+
+In order to include the spin-orbit interaction into the scalar-relativistic formalism of a normconserving, non-local pseudopotential, scalar-relativistic in contrast to fully relativistic pseudopotential wavefunctions are needed as a basis for the projectors of the spin-orbit potential.
+The transformation
+\begin{equation}
+L\cdot S=L_xS_x+L_yS_y+L_zS_z
+\end{equation}
+using the ladder operators
+\begin{equation}
+L_\pm=L_x\pm iL_y \text{ and } S_\pm=S_x\pm iS_y
+\text{ ,}
+\end{equation}
+with properties
+\begin{eqnarray}
+L_+S_- & = & (L_x+iL_y)(S_x-iS_y)=L_xS_x-iL_xS_y+iL_yS_x+L_yS_y \\
+L_-S_+ & = & (L_x-iL_y)(S_x+iS_y)=L_xS_x+iL_xS_y-iL_yS_x+L_yS_y 
+\end{eqnarray}
+resulting in
+\begin{equation}
+L_+S_-+L_-S_+=2(L_xS_x+L_yS_y)
+\text{ ,}
+\end{equation}
+reads
+\begin{equation}
+L\cdot S=\frac{1}{2}(L_+S_-+L_-S_+)+L_zS_z
+\text{ .}
+\end{equation}
+The contributions of this operator act differently on $\ket{l,m}$ and --- in fact --- depend on the respectively considered spinor component, which is incorporated by $\ket{l,m,\pm}$.
+\begin{enumerate}
+\item \underline{$L_+S_-$}:
+      Updates spin down component and only acts on spin up component
+\begin{equation}
+L_+S_-\ket{l,m,+}=L_+\ket{l,m}S_-\ket{+}=
+\sqrt{(l-m)(l+m+1)}\hbar\ket{l,m+1}\hbar\ket{-}
+\end{equation}
+\item \underline{$L_-S_+$}:
+      Updates spin up component and only acts on spin down component
+\begin{equation}
+L_+S_-\ket{l,m,-}=L_+\ket{l,m}S_-\ket{+}=
+\sqrt{(l-m)(l+m+1)}\hbar\ket{l,m+1}\hbar\ket{-}
+\end{equation}
+\item \underline{$L_zS_z$}: Acts on both and updates both spinor components
+\begin{equation}
+L_zS_z\ket{l,m,\pm}=L_z\ket{l,m}S_z\ket{\pm}=
+\pm\frac{1}{2}m\hbar^2\ket{l,m,\pm}
+\end{equation}
+\end{enumerate}
+
+\subsubsection{Excursus: Real space representation within an iterative treatment}
+
+In the following, the spin-orbit part is evaluated in real space.
+Since spin is treated in another subspace, it can be treated separately.
+The matrix elements of the orbital angular momentum part of the potential in KB form read
+\begin{equation}
+\sum_{lm}
+\bra{\vec{r}'}(r\times p)\ket{\chi_{lm}}E^{\text{SO,KB}}_l
+\braket{\chi_{lm}}{\vec{r}''}
+\text{ .}
+\end{equation}
+With
+\begin{eqnarray}
+\bra{\vec{r}'}r\ket{\chi_{lm}} & = & \vec{r}'\braket{\vec{r}'}{\chi_{lm}}\\
+\bra{\vec{r}'}p\ket{\chi_{lm}} & = & -i\hbar\nabla_{\vec{r}'}
+\braket{\vec{r}'}{\chi_{lm}}
+\end{eqnarray}
+we get
+\begin{equation}
+-i\hbar\sum_{lm}(\vec{r}'\times \nabla_{\vec{r}'})\braket{\vec{r}'}{\chi_{lm}}
+E^{\text{SO,KB}}_l\braket{\chi_{lm}}{\vec{r}''}
+\text{ .}
+\label{eq:solid:so_me}
+\end{equation}
+To further evaluate this expression, the KB projectors
+\begin{equation}
+\ket{\chi_{lm}}=\frac{\ket{\delta V_l^{\text{SO}}\Phi_{lm}}}
+{\braket{\delta V_l^{\text{SO}}\Phi_{lm}}
+        {\Phi_{lm}\delta V_l^{\text{SO}}}^{1/2}}
+\end{equation}
+must be known in real space (with respect to $\vec{r}'$).
+\begin{equation}
+\braket{\vec{r}'}{\chi_{lm}}=
+\frac{\braket{\vec{r}'}{\delta V_l^{\text{SO}}\Phi_{lm}}}{
+\braket{\delta V_l^{\text{SO}}\Phi_{lm}}{\Phi_{lm}\delta V_l^{\text{SO}}}
+^{1/2}}
+\end{equation}
+and
+\begin{equation}
+\braket{\vec{r}'}{\delta V_l^{\text{SO}}\Phi_{lm}}=
+\delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}Y_{lm}(\Omega_{\vec{r}'})
+\text{ .}
+\label{eq:solid:so_r1}
+\end{equation}
+In this expression, only the spherical harmonics are complex functions.
+Thus, the complex conjugate with respect to $\vec{r}''$ is given by
+\begin{equation}
+\braket{\Phi_{lm}\delta V_l^{\text{SO}}}{\vec{r}''}=
+\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}Y^*_{lm}(\Omega_{\vec{r}''})
+\text{ .}
+\label{eq:solid:so_r2}
+\end{equation}
+Using the orthonormality property 
+\begin{equation}
+\int Y^*_{l'm'}(\Omega_r)Y_{lm}(\Omega_r) d\Omega_r = \delta_{ll'}\delta_{mm'}
+\label{eq:solid:y_ortho}
+\end{equation}
+of the spherical harmonics, the squared norm of the $\chi_{lm}$ reduces to
+\begin{eqnarray}
+\braket{\delta V_l^{\text{SO}}\Phi_{lm}}{\Phi_{lm}\delta V_l^{\text{SO}}} &=&
+\int \braket{\delta V_l^{\text{SO}}\Phi_{lm}}{\vec{r}'}
+\braket{\vec{r}'}{\Phi_{lm}\delta V_l^{\text{SO}}} d\vec{r}'\\
+&=&\int 
+{\delta V_l^{\text{SO}}}^2(r')\frac{u_l^2(r')}{r'^2}Y^*_{lm}(\Omega_{\vec{r}'})
+Y_{lm}(\Omega_{\vec{r}'})
+r'^2 dr' d\Omega_{\vec{r}'} \\
+&=&\int_{r'}
+{\delta V_l^{\text{SO}}}^2(r') u_l^2(r') dr'
+\int_{\Omega_{\vec{r}'}}Y^*_{lm}(\Omega_{\vec{r}'})Y_{lm}(\Omega_{\vec{r}'}) d\Omega_{\vec{r}'}\\
+&=&\int_{r'}{\delta V_l^{\text{SO}}}^2(r') u_l^2(r') dr' \text{ .}\\
+&=&\braket{\delta V_l^{\text{SO}}u_l}{u_l\delta V_l^{\text{SO}}}
+\end{eqnarray}
+To obtain a final expression for the matrix elements \eqref{eq:solid:so_me}, the sum of the products of \eqref{eq:solid:so_r1} and \eqref{eq:solid:so_r2} must be further evaluated.
+\begin{eqnarray}
+\sum_{lm}
+\braket{\vec{r}'}{\delta V_l^{\text{SO}}\Phi_{lm}}
+\braket{\Phi_{lm}\delta V_l^{\text{SO}}}{\vec{r}''}&=&\sum_{lm}
+\delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}Y_{lm}(\Omega_{\vec{r}'})
+\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}Y^*_{lm}(\Omega_{\vec{r}''})\nonumber\\
+&=&\sum_l
+\delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
+\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}\sum_m
+Y^*_{lm}(\Omega_{\vec{r}''})Y_{lm}(\Omega_{\vec{r}'})\nonumber\\
+&=&\sum_l
+\delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
+\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}
+P_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)\frac{2l+1}{4\pi}\nonumber\\
+\end{eqnarray}
+due to the vector addition theorem
+\begin{equation}
+P_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)=
+\frac{4\pi}{2l+1}\sum_mY^*_{lm}(\Omega_{\vec{r}''})Y_{lm}(\Omega_{\vec{r}'})
+\text{ .}
+\end{equation}
+In total, the matrix elements of the SO potential can be calculated by
+\begin{eqnarray}
+&&-i\hbar\sum_{lm}(\vec{r}'\times \nabla_{\vec{r}'})\braket{\vec{r}'}{\chi_{lm}}
+E^{\text{SO,KB}}_l\braket{\chi_{lm}}{\vec{r}''}=\nonumber\\
+&=&-i\hbar\sum_l(\vec{r}'\times \nabla_{\vec{r}'})
+\delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
+P_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)\cdot
+\frac{E^{\text{SO,KB}}_l\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}}
+             {\int_{r}{\delta V_l^{\text{SO}}}^2(r) u_l^2(r) dr} \cdot
+\frac{2l+1}{4\pi}\nonumber\\
+&=&
+-i\hbar\sum_l
+\delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
+P'_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)\cdot
+\left(\frac{\vec{r}'\times\vec{r}''}{r'r''}\right)\cdot
+\frac{E^{\text{SO,KB}}_l\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}}
+       {\int_{r}{\delta V_l^{\text{SO}}}^2(r) u_l^2(r) dr} \,
+\frac{2l+1}{4\pi}\text{ ,}\nonumber\\
+\label{eq:solid:so_fin}
+\end{eqnarray}
+since derivatives of functions that depend on the absolute value of $\vec{r}'$ do not contribute due to the cross product as is illustrated below (equations \eqref{eq:solid:rxp1} and \eqref{eq:solid:rxp2}).
+\begin{eqnarray}
+\left(\vec{r}\times\nabla_{\vec{r}}\right)f(r)&=&
+\left(\begin{array}{l}
+r_y\frac{\partial}{\partial r_z}f(r)-r_z\frac{\partial}{\partial r_y}f(r)\\
+r_z\frac{\partial}{\partial r_x}f(r)-r_x\frac{\partial}{\partial r_z}f(r)\\
+r_x\frac{\partial}{\partial r_y}f(r)-r_y\frac{\partial}{\partial r_x}f(r)
+\end{array}\right)
+\label{eq:solid:rxp1}
+\end{eqnarray}
+\begin{eqnarray}
+r_i\frac{\partial}{\partial r_j}f(r)-r_j\frac{\partial}{\partial r_i}f(r)&=&
+r_if'(r)\frac{\partial}{\partial r_j}(r_x^2+r_y^2+r_z^2)^{1/2}-
+r_jf'(r)\frac{\partial}{\partial r_i}(r_x^2+r_y^2+r_z^2)^{1/2}\nonumber\\
+&=&
+r_if'(r)\frac{1}{2r}2r_j-r_jf'(r)\frac{1}{2r}2r_i=0
+\label{eq:solid:rxp2}
+\end{eqnarray}
+
+If these projectors are considered to be centered around atom positions $\vec{\tau}_{\alpha n}$ of atoms $n$ of species $\alpha$, the variable $\vec{r}'$ in the previous equations is changed to $\vec{r}'_{\alpha n}=\vec{r}'-\vec{\tau}_{\alpha n}$, which implies
+\begin{eqnarray}
+r'&\rightarrow&r_{\alpha n}=|\vec{r}'-\vec{\tau}_{\alpha n}|\\
+\Omega_{\vec{r}'}&\rightarrow&\Omega_{\vec{r'}-\vec{\tau}_{\alpha n}}\\
+\delta V_l(r')&\rightarrow&\delta V_l(|\vec{r}'-\vec{\tau}_{\alpha n}|)\\
+u_l(r')&\rightarrow&u_l(|\vec{r}'-\vec{\tau}_{\alpha n}|)\\
+Y_{lm}(\Omega_{\vec{r}'})&\rightarrow&
+Y_{lm}(\Omega_{\vec{r}'-\vec{\tau}_{\alpha n}})
+\text{ .}
+\end{eqnarray}
+Within an iterative treatment on a real space grid consisting of $n_{\text{g}}$ grid points, the sum
+\begin{equation}
+\sum_{\vec{r}''_{\alpha n}}
+\sum_{lm}-i\hbar(\vec{r}'_{\alpha n}\times \nabla_{\vec{r}'_{\alpha n}})
+\braket{\vec{r}'_{\alpha n}}{\chi^{\text{SO}}_{lm}}
+E^{\text{SO,KB}}_l\braket{\chi^{\text{SO}}_{lm}}{\vec{r}''_{\alpha n}}
+\braket{\vec{r}''_{\alpha n}}{\Psi}
+\qquad\forall\,\bra{\vec{r}'_{\alpha n}}
+\end{equation}
+to obtain all elements $\bra{\vec{r}'_{\alpha n}}$, involves $n_{\text{g}}^2$ evaluations of equation~\eqref{eq:solid:so_fin} for eeach atom, if the projectors are short-ranged, i.e.\  $\delta V_l=0$ outside a certain cut-off radius.
+Thus, this method scales linearly with the number of atoms.
+
+The $E_l^{\text{SO,KB}}$ are given by
+\begin{equation}
+E_l^{\text{SO,KB}}=
+\frac{\braket{\delta V_lu_l}{u_l\delta V_l}}
+     {\bra{u_l}\delta V_l\ket{u_l}}=
+\frac{\int_{r}\delta V^2_l(r)u^2_l(r)}r^2dr
+     {\int_{r'}\int_{r''}\braket{u_l}{r'}\bra{r'}\delta V_l
+\ket{r''}\braket{r''}{u_l}}=
 \end{equation}
-By switching the indices ...