physics_compact/solid.tex

   1 \part{Theory of the solid state}
   2
   3 \chapter{Atomic structure}
   4
   5 \chapter{Reciprocal lattice}
   6
   7 Example of primitive cell ...
   8
   9 \chapter{Electronic structure}
  10
  11 \section{Noninteracting electrons}
  12
  13 \subsection{Bloch's theorem}
  14
  15 \section{Nearly free and tightly bound electrons}
  16
  17 \subsection{Tight binding model}
  18
  19 \section{Interacting electrons}
  20
  21 \subsection{Density functional theory}
  22
  23 \subsubsection{Hohenberg-Kohn theorem}
  24
  25 The Hamiltonian of a many-electron problem has the form
  26 \begin{equation}
  27 H=T+V+U\text{ ,}
  28 \end{equation}
  29 where
  30 \begin{eqnarray}
  31 T & = & \langle\Psi|\sum_{i=1}^N\frac{-\hbar^2}{2m}\nabla_i^2|\Psi\rangle\\
  32   & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} d\vec{r}' \,
  33         \langle \Psi | \vec{r} \rangle \langle \vec{r} |
  34         \nabla_i^2
  35         | \vec{r}' \rangle \langle \vec{r}' | \Psi \rangle\\
  36   & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} d\vec{r}' \,
  37         \langle \Psi | \vec{r} \rangle \nabla_{\vec{r}_i}
  38         \langle \vec{r} | \vec{r}' \rangle
  39         \nabla_{\vec{r}'_i} \langle \vec{r}' | \Psi \rangle\\
  40   & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} d\vec{r}' \,
  41         \nabla_{\vec{r}_i} \langle \Psi | \vec{r} \rangle
  42         \delta_{\vec{r}\vec{r}'}
  43         \nabla_{\vec{r}'_i} \langle \vec{r}' | \Psi \rangle\\
  44   & = & \frac{-\hbar^2}{2m} \sum_{i=1}^N \int d\vec{r} \,
  45         \nabla_{\vec{r}_i} \Psi^*(\vec{r}) \nabla_{\vec{r}_i} \Psi(\vec{r})
  46         \text{ ,} \\
  47 V & = & \int V(\vec{r})\Psi^*(\vec{r})\Psi(\vec{r})d\vec{r} \text{ ,} \\
  48 U & = & \frac{1}{2}\int\frac{1}{\left|\vec{r}-\vec{r}'\right|}
  49         \Psi^*(\vec{r})\Psi^*(\vec{r}')\Psi(\vec{r}')\Psi(\vec{r})
  50         d\vec{r}d\vec{r}'
  51 \end{eqnarray}
  52 represent the kinetic energy, the energy due to the external potential and the energy due to the mutual Coulomb repulsion.
  53
  54 \begin{remark}
  55 As can be seen from the above, two many-electron systems can only differ in the external potential and the number of electrons.
  56 The number of electrons is determined by the electron density.
  57 \begin{equation}
  58 N=\int n(\vec{r})d\vec{r}
  59 \end{equation}
  60 Now, if the external potential is additionally determined by the electron density, the density completely determines the many-body problem.
  61 \end{remark}
  62
  63 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})$.
  64 \begin{equation}
  65 n_0(\vec{r})=\int \Psi_0^*(\vec{r},\vec{r}_2,\vec{r}_3,\ldots,\vec{r}_N)
  66                   \Psi_0(\vec{r},\vec{r}_2,\vec{r}_3,\ldots,\vec{r}_N)
  67              d\vec{r}_2d\vec{r}_3\ldots d\vec{r}_N
  68 \end{equation}
  69 In 1964, Hohenberg and Kohn showed the opposite and far less obvious result \cite{hohenberg64}.
  70
  71 \begin{theorem}[Hohenberg / Kohn]
  72 For a nondegenerate ground state, the ground-state charge density uniquely determines the external potential in which the electrons reside.
  73 \end{theorem}
  74
  75 \begin{proof}
  76 The proof presented by Hohenberg and Kohn proceeds by {\em reductio ad absurdum}.
  77 Suppose two potentials $V_1$ and $V_2$ exist, which yield the same electron density $n(\vec{r})$.
  78 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$.
  79 Then, due to the variational principle (see \ref{sec:var_meth}), one can write
  80 \begin{equation}
  81 E_1=\langle \Psi_1 | H_1 | \Psi_1 \rangle <
  82 \langle \Psi_2 | H_1 | \Psi_2 \rangle \text{ .}
  83 \label{subsub:hk01}
  84 \end{equation}
  85 Expressing $H_1$ by $H_2+H_1-H_2$, the last part of \eqref{subsub:hk01} can be rewritten:
  86 \begin{equation}
  87 \langle \Psi_2 | H_1 | \Psi_2 \rangle =
  88 \langle \Psi_2 | H_2 | \Psi_2 \rangle +
  89 \langle \Psi_2 | H_1 -H_2 | \Psi_2 \rangle
  90 \end{equation}
  91 The two Hamiltonians, which describe the same number of electrons, differ only in the potential
  92 \begin{equation}
  93 H_1-H_2=V_1(\vec{r})-V_2(\vec{r})
  94 \end{equation}
  95 and, thus
  96 \begin{equation}
  97 E_1<E2+\int n(\vec{r}) \left( V_1(\vec{r})-V_2(\vec{r}) \right) d\vec{r}
  98 \text{ .}
  99 \label{subsub:hk02}
 100 \end{equation}
 101 By switching the indices of \eqref{subsub:hk02} and adding the resulting equation to \eqref{subsub:hk02}, the contradiction
 102 \begin{equation}
 103 E_1 + E_2 < E_2 + E_1 +
 104 \underbrace{
 105 \int n(\vec{r}) \left( V_1(\vec{r})-V_2(\vec{r}) \right) d\vec{r} +
 106 \int n(\vec{r}) \left( V_2(\vec{r})-V_1(\vec{r}) \right) d\vec{r}
 107 }_{=0}
 108 \end{equation}
 109 is revealed, which proofs the Hohenberg Kohn theorem.% \qed
 110 \end{proof}
 111
 112 \section{Electron-ion interaction}
 113
 114 \subsection{Pseudopotential theory}
 115
 116 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.
 117
 118 \subsubsection{Orthogonalized planewave method}
 119
 120 Due to the orthogonality of the core and valence wavefunctions, the latter exhibit strong oscillations within the core region of the atom.
 121 This requires a large amount of planewaves $\ket{k}$ to adequatley model the valence electrons.
 122
 123 In a very general approach, the orthogonalized planewave (OPW) method introduces a new basis set
 124 \begin{equation}
 125 \ket{k}_{\text{OPW}} = \ket{k} - \sum_t \ket{t}\bra{t}k\rangle \text{ ,}
 126 \end{equation}
 127 with $\ket{t}$ being the eigenstates of the core electrons.
 128 The new basis is orthogonal to the core states $\ket{t}$.
 129 \begin{equation}
 130 \braket{t}{k}_{\text{OPW}} =
 131 \braket{t}{k} - \sum_{t'} \braket{t}{t'}\braket{t'}{k} =
 132 \braket{t}{k} - \braket{t}{k}=0
 133 \end{equation}
 134 The number of planewaves required for reasonably converged electronic structure calculations is tremendously reduced by utilizing the OPW basis set.
 135
 136 \subsubsection{Pseudopotential method}
 137
 138 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.
 139
 140 Let $\ket{\Psi_\text{V}}$ be the wavefunction of a valence electron with the Schr\"odinger equation
 141 \begin{equation}
 142 H \ket{\Psi_\text{V}} = \left(\frac{1}{2m}p^2+V\right)\ket{\Psi_\text{V}}=
 143 E\ket{\Psi_\text{V}} \text{ .}
 144 \end{equation}
 145 \ldots projection operatore $P_\text{C}$ \ldots
 146
 147 \subsubsection{Semilocal form of the pseudopotential}
 148
 149 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
 150 \begin{equation}
 151 V=\sum_{l,m}\ket{lm}V_l(\vec{r})\bra{lm} \text{ .}
 152 \end{equation}
 153 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.
 154
 155 Problem of semilocal potantials become valid once matrix elements need to be computed.
 156 Integral with respect to the radial component needs to be evaluated for each planewave combination, i.e.\  $N(N-1)/2$ integrals.
 157 \begin{equation}
 158 \bra{k+G}V\ket{k+G'} = \ldots
 159 \end{equation}
 160
 161 A local potential can always be separated from the potential \ldots
 162 \begin{equation}
 163 V=\ldots=V_{\text{local}}(\vec{r})+\ldots
 164 \end{equation}
 165
 166 \subsubsection{Norm conserving pseudopotentials}
 167
 168 HSC potential \ldots
 169
 170 \subsubsection{Fully separable form of the pseudopotential}
 171
 172 KB transformation \ldots
 173
 174 \subsection{Spin-orbit interaction}
 175
 176 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.
 177 This is advantageous since \ldots
 178 With the solutions of the all-electron Dirac equations, the new pseudopotential reads
 179 \begin{equation}
 180 V(\vec{r})=\sum_{l,m}\left[
 181 \ket{l+\frac{1}{2},m+{\frac{1}{2}}}V_{l,l+\frac{1}{2}}(\vec{r})
 182 \bra{l+\frac{1}{2},m+{\frac{1}{2}}} +
 183 \ket{l-\frac{1}{2},m-{\frac{1}{2}}}V_{l,l-\frac{1}{2}}(\vec{r})
 184 \bra{l-\frac{1}{2},m-{\frac{1}{2}}}
 185 \right] \text{ .}
 186 \label{eq:solid:so_bs1}
 187 \end{equation}
 188 By defining an averaged potential weighted by the different $j$ degeneracies of the $\ket{l\pm\frac{1}{2}}$ states
 189 \begin{equation}
 190 \bar{V}_l(r)=\frac{1}{2l+1}\left(
 191 l V_{l,l-\frac{1}{2}}(\vec{r})+(l+1)V_{l,l+\frac{1}{2}}(\vec{r})\right)
 192 \end{equation}
 193 and a potential describing the difference in the potential with respect to the spin
 194 \begin{equation}
 195 V^{\text{SO}}_l(\vec{r})=\frac{2}{2l+1}\left(
 196 V_{l,l+\frac{1}{2}}(\vec{r})-V_{l,l-\frac{1}{2}}(\vec{r})\right)
 197 \end{equation}
 198 the total potential can be expressed as
 199 \begin{equation}
 200 V(\vec{r})=\sum_l
 201 \ket{l,m}\left[\bar{V}_l(\vec{r})+V^{\text{SO}}_l(\vec{r})LS\right]\bra{l,m}
 202 \text{ ,}
 203 \label{eq:solid:so_bs2}
 204 \end{equation}
 205 where the first term correpsonds to the mass velocity and Darwin relativistic corrections and the latter is associated with the spin-orbit (SO) coupling.
 206 \begin{proof}
 207 This can be shown by rewriting the $LS$ operator
 208 \begin{equation}
 209 J=L+S \Leftrightarrow J^2=L^2+S^2+2LS \Leftrightarrow
 210 LS=\frac{1}{2}\left(J^2-L^2-S^2\right)
 211 \end{equation}
 212 and corresponding eigenvalue
 213 \begin{eqnarray}
 214 j(j+1)-l(l+1)-s(s+1)&=&
 215 (l\pm\frac{1}{2})(l\pm\frac{1}{2}+1)-l^2-l-\frac{3}{4} \nonumber\\
 216 &=&
 217 l^2\pm\frac{l}{2}+l\pm\frac{l}{2}+\frac{1}{4}\pm\frac{1}{2}-l^2-l-\frac{3}{4}
 218 \nonumber\\
 219 &=&\pm(l+\frac{1}{2})-\frac{1}{2}=\left\{\begin{array}{rl}
 220 l & \text{for } j=l+\frac{1}{2}\\
 221 -(l+1) & \text{for } j=l-\frac{1}{2}
 222 \end{array}\right.
 223 \text{ ,}
 224 \end{eqnarray}
 225 which, if used in equation~\eqref{eq:solid:so_bs2}, gives the same (diagonal) matrix elements
 226 \begin{eqnarray}
 227 \bra{l\pm\frac{1}{2},m\pm\frac{1}{2}}V(\vec{r})
 228 \ket{l\pm\frac{1}{2},m\pm\frac{1}{2}}&=&
 229 \bar{V}_l(\vec{r})+V^{\text{SO}}_l(\vec{r})
 230 \frac{1}{2}\left(l(l+1)-j(j+1)-\frac{3}{4}\right) \nonumber\\
 231 &=&\bar{V}_l(\vec{r})+\frac{1}{2}V^{\text{SO}}_l(\vec{r})
 232 \cdot\left\{\begin{array}{cl}
 233 l & \text{for } j=l+\frac{1}{2}\\
 234 -(l+1) & \text{for } j=l-\frac{1}{2}
 235 \end{array}\right. \nonumber\\
 236 &=&\frac{1}{2l+1}\left(lV_{l,l-\frac{1}{2}}(\vec{r})+
 237                        (l+1)V_{l,l+\frac{1}{2}}(\vec{r})\right)+\nonumber\\
 238 &&\frac{1}{2l+1}
 239 \left(V_{l,l+\frac{1}{2}}(\vec{r})-V_{l,l-\frac{1}{2}}(\vec{r})\right)
 240 \cdot\left\{\begin{array}{c}
 241 l \\
 242 -(l+1)
 243 \end{array}\right. \nonumber\\
 244 &=&\left\{\begin{array}{cl}
 245 V_{l,l+\frac{1}{2}}(\vec{r}) & \text{for } j=l+\frac{1}{2}\\
 246 V_{l,l-\frac{1}{2}}(\vec{r}) & \text{for } j=l-\frac{1}{2}
 247 \end{array}\right.
 248 \end{eqnarray}
 249 as expected and --- in fact --- obtained from equation~\eqref{eq:solid:so_bs1}.
 250 \end{proof}
 251
 252 \subsubsection{Scalar relativistic basis}
 253
 254 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.
 255 The transformation
 256 \begin{equation}
 257 L\cdot S=L_xS_x+L_yS_y+L_zS_z
 258 \end{equation}
 259 using the ladder operators
 260 \begin{equation}
 261 L_\pm=L_x\pm iL_y \text{ and } S_\pm=S_x\pm iS_y
 262 \text{ ,}
 263 \end{equation}
 264 with properties
 265 \begin{eqnarray}
 266 L_+S_- & = & (L_x+iL_y)(S_x-iS_y)=L_xS_x-iL_xS_y+iL_yS_x+L_yS_y \\
 267 L_-S_+ & = & (L_x-iL_y)(S_x+iS_y)=L_xS_x+iL_xS_y-iL_yS_x+L_yS_y
 268 \end{eqnarray}
 269 resulting in
 270 \begin{equation}
 271 L_+S_-+L_-S_+=2(L_xS_x+L_yS_y)
 272 \text{ ,}
 273 \end{equation}
 274 reads
 275 \begin{equation}
 276 L\cdot S=\frac{1}{2}(L_+S_-+L_-S_+)+L_zS_z
 277 \text{ .}
 278 \end{equation}
 279 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}$.
 280 \begin{enumerate}
 281 \item \underline{$L_+S_-$}:
 282       Updates spin down component and only acts on spin up component
 283 \begin{equation}
 284 L_+S_-\ket{l,m,+}=L_+\ket{l,m}S_-\ket{+}=
 285 \sqrt{(l-m)(l+m+1)}\hbar\ket{l,m+1}\hbar\ket{-}
 286 \end{equation}
 287       Moreover, this part only acts on magnetic quantum numbers
 288       $m=-l,\ldots,l-1$ and updates quantum numbers $m=-l+1,\ldots,l$.
 289 \item \underline{$L_-S_+$}:
 290       Updates spin up component and only acts on spin down component
 291 \begin{equation}
 292 L_-S_+\ket{l,m,-}=L_+\ket{l,m}S_+\ket{-}=
 293 \sqrt{(l+m)(l-m+1)}\hbar\ket{l,m-1}\hbar\ket{+}
 294 \end{equation}
 295       Moreover, this part only acts on magnetic quantum numbers
 296       $m=-l+1,\ldots,l$ and updates quantum numbers $m=-l,\ldots,l-1$.
 297 \item \underline{$L_zS_z$}: Acts on both and updates both spinor components
 298 \begin{equation}
 299 L_zS_z\ket{l,m,\pm}=L_z\ket{l,m}S_z\ket{\pm}=
 300 \pm\frac{1}{2}m\hbar^2\ket{l,m,\pm}
 301 \end{equation}
 302       It acts on all magnetic quantum numbers and updates all of them.
 303 \end{enumerate}
 304 Please note that the $\ket{l,m,\pm}$ are not eigenfunctions of the two combinations of ladder operators, i.e.\  the $\ket{l,m,\pm}$ do not diagonalize the spin-orbit part of the Hamiltonian.
 305
 306 These equations can be simplified to read
 307 \begin{eqnarray}
 308 \ldots
 309 \text{ .}
 310 \end{eqnarray}
 311
 312 \subsubsection{A different basis set}
 313
 314 The above basis is composed of eigenfunctions
 315 \begin{equation}
 316 \ket{l,m} \text{, } \ket{\pm} \text{ of operators }
 317 L^2\text{, } L_z \text{ and } S_z
 318 \text{.}
 319 \end{equation}
 320 These eigenfunctions diagonalize the scalar relativistic Hamiltonian.
 321 Introducing spin-orbit interaction, however, it is a good idea to chose eigenfunctions that diagonalize the perturbation
 322 \begin{equation}
 323 L\cdot S=\frac{1}{2}(J^2-L^2-S^2)
 324 \text{ ,}
 325 \end{equation}
 326 i.e.\  simultaneous eigenfunctions of $J^2$, $L^2$ and $S^2$.
 327
 328 \subsubsection{Excursus: Real space representation within an iterative treatment}
 329
 330 In the following, the spin-orbit part is evaluated in real space.
 331 Since spin is treated in another subspace, it can be treated separately.
 332 The matrix elements of the orbital angular momentum part of the potential in KB form read
 333 \begin{equation}
 334 \sum_{lm}
 335 \bra{\vec{r}'}(r\times p)\ket{\chi_{lm}}E^{\text{SO,KB}}_l
 336 \braket{\chi_{lm}}{\vec{r}''}
 337 \text{ .}
 338 \end{equation}
 339 With
 340 \begin{eqnarray}
 341 \bra{\vec{r}'}r\ket{\chi_{lm}} & = & \vec{r}'\braket{\vec{r}'}{\chi_{lm}}\\
 342 \bra{\vec{r}'}p\ket{\chi_{lm}} & = & -i\hbar\nabla_{\vec{r}'}
 343 \braket{\vec{r}'}{\chi_{lm}}
 344 \end{eqnarray}
 345 we get
 346 \begin{equation}
 347 -i\hbar\sum_{lm}(\vec{r}'\times \nabla_{\vec{r}'})\braket{\vec{r}'}{\chi_{lm}}
 348 E^{\text{SO,KB}}_l\braket{\chi_{lm}}{\vec{r}''}
 349 \text{ .}
 350 \label{eq:solid:so_me}
 351 \end{equation}
 352 To further evaluate this expression, the KB projectors
 353 \begin{equation}
 354 \ket{\chi_{lm}}=\frac{\ket{\delta V_l^{\text{SO}}\Phi_{lm}}}
 355 {\braket{\delta V_l^{\text{SO}}\Phi_{lm}}
 356         {\Phi_{lm}\delta V_l^{\text{SO}}}^{1/2}}
 357 \end{equation}
 358 must be known in real space (with respect to $\vec{r}'$).
 359 \begin{equation}
 360 \braket{\vec{r}'}{\chi_{lm}}=
 361 \frac{\braket{\vec{r}'}{\delta V_l^{\text{SO}}\Phi_{lm}}}{
 362 \braket{\delta V_l^{\text{SO}}\Phi_{lm}}{\Phi_{lm}\delta V_l^{\text{SO}}}
 363 ^{1/2}}
 364 \end{equation}
 365 and
 366 \begin{equation}
 367 \braket{\vec{r}'}{\delta V_l^{\text{SO}}\Phi_{lm}}=
 368 \delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}Y_{lm}(\Omega_{\vec{r}'})
 369 \text{ .}
 370 \label{eq:solid:so_r1}
 371 \end{equation}
 372 In this expression, only the spherical harmonics are complex functions.
 373 Thus, the complex conjugate with respect to $\vec{r}''$ is given by
 374 \begin{equation}
 375 \braket{\Phi_{lm}\delta V_l^{\text{SO}}}{\vec{r}''}=
 376 \delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}Y^*_{lm}(\Omega_{\vec{r}''})
 377 \text{ .}
 378 \label{eq:solid:so_r2}
 379 \end{equation}
 380 Using the orthonormality property
 381 \begin{equation}
 382 \int Y^*_{l'm'}(\Omega_r)Y_{lm}(\Omega_r) d\Omega_r = \delta_{ll'}\delta_{mm'}
 383 \label{eq:solid:y_ortho}
 384 \end{equation}
 385 of the spherical harmonics, the squared norm of the $\chi_{lm}$ reduces to
 386 \begin{eqnarray}
 387 \braket{\delta V_l^{\text{SO}}\Phi_{lm}}{\Phi_{lm}\delta V_l^{\text{SO}}} &=&
 388 \int \braket{\delta V_l^{\text{SO}}\Phi_{lm}}{\vec{r}'}
 389 \braket{\vec{r}'}{\Phi_{lm}\delta V_l^{\text{SO}}} d\vec{r}'\\
 390 &=&\int
 391 {\delta V_l^{\text{SO}}}^2(r')\frac{u_l^2(r')}{r'^2}Y^*_{lm}(\Omega_{\vec{r}'})
 392 Y_{lm}(\Omega_{\vec{r}'})
 393 r'^2 dr' d\Omega_{\vec{r}'} \\
 394 &=&\int_{r'}
 395 {\delta V_l^{\text{SO}}}^2(r') u_l^2(r') dr'
 396 \int_{\Omega_{\vec{r}'}}Y^*_{lm}(\Omega_{\vec{r}'})Y_{lm}(\Omega_{\vec{r}'}) d\Omega_{\vec{r}'}\\
 397 &=&\int_{r'}{\delta V_l^{\text{SO}}}^2(r') u_l^2(r') dr' \text{ .}\\
 398 &=&\braket{\delta V_l^{\text{SO}}u_l}{u_l\delta V_l^{\text{SO}}}
 399 \end{eqnarray}
 400 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.
 401 \begin{eqnarray}
 402 \sum_{lm}
 403 \braket{\vec{r}'}{\delta V_l^{\text{SO}}\Phi_{lm}}
 404 \braket{\Phi_{lm}\delta V_l^{\text{SO}}}{\vec{r}''}&=&\sum_{lm}
 405 \delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}Y_{lm}(\Omega_{\vec{r}'})
 406 \delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}Y^*_{lm}(\Omega_{\vec{r}''})\nonumber\\
 407 &=&\sum_l
 408 \delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
 409 \delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}\sum_m
 410 Y^*_{lm}(\Omega_{\vec{r}''})Y_{lm}(\Omega_{\vec{r}'})\nonumber\\
 411 &=&\sum_l
 412 \delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
 413 \delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}
 414 P_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)\frac{2l+1}{4\pi}\nonumber\\
 415 \end{eqnarray}
 416 due to the vector addition theorem
 417 \begin{equation}
 418 P_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)=
 419 \frac{4\pi}{2l+1}\sum_mY^*_{lm}(\Omega_{\vec{r}''})Y_{lm}(\Omega_{\vec{r}'})
 420 \text{ .}
 421 \end{equation}
 422 In total, the matrix elements of the SO potential can be calculated by
 423 \begin{eqnarray}
 424 &&-i\hbar\sum_{lm}(\vec{r}'\times \nabla_{\vec{r}'})\braket{\vec{r}'}{\chi_{lm}}
 425 E^{\text{SO,KB}}_l\braket{\chi_{lm}}{\vec{r}''}=\nonumber\\
 426 &=&-i\hbar\sum_l(\vec{r}'\times \nabla_{\vec{r}'})
 427 \delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
 428 P_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)\cdot
 429 \frac{E^{\text{SO,KB}}_l\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}}
 430              {\int_{r}{\delta V_l^{\text{SO}}}^2(r) u_l^2(r) dr} \cdot
 431 \frac{2l+1}{4\pi}\nonumber\\
 432 &=&
 433 -i\hbar\sum_l
 434 \delta V_l^{\text{SO}}(r')\frac{u_l(r')}{r'}
 435 P'_l\left(\frac{\vec{r}'\vec{r}''}{r'r''}\right)\cdot
 436 \left(\frac{\vec{r}'\times\vec{r}''}{r'r''}\right)\cdot
 437 \frac{E^{\text{SO,KB}}_l\delta V_l^{\text{SO}}(r'')\frac{u_l(r'')}{r''}}
 438        {\int_{r}{\delta V_l^{\text{SO}}}^2(r) u_l^2(r) dr} \,
 439 \frac{2l+1}{4\pi}\text{ ,}\nonumber\\
 440 \label{eq:solid:so_fin}
 441 \end{eqnarray}
 442 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}).
 443 \begin{eqnarray}
 444 \left(\vec{r}\times\nabla_{\vec{r}}\right)f(r)&=&
 445 \left(\begin{array}{l}
 446 r_y\frac{\partial}{\partial r_z}f(r)-r_z\frac{\partial}{\partial r_y}f(r)\\
 447 r_z\frac{\partial}{\partial r_x}f(r)-r_x\frac{\partial}{\partial r_z}f(r)\\
 448 r_x\frac{\partial}{\partial r_y}f(r)-r_y\frac{\partial}{\partial r_x}f(r)
 449 \end{array}\right)
 450 \label{eq:solid:rxp1}
 451 \end{eqnarray}
 452 \begin{eqnarray}
 453 r_i\frac{\partial}{\partial r_j}f(r)-r_j\frac{\partial}{\partial r_i}f(r)&=&
 454 r_if'(r)\frac{\partial}{\partial r_j}(r_x^2+r_y^2+r_z^2)^{1/2}-
 455 r_jf'(r)\frac{\partial}{\partial r_i}(r_x^2+r_y^2+r_z^2)^{1/2}\nonumber\\
 456 &=&
 457 r_if'(r)\frac{1}{2r}2r_j-r_jf'(r)\frac{1}{2r}2r_i=0
 458 \label{eq:solid:rxp2}
 459 \end{eqnarray}
 460
 461 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
 462 \begin{eqnarray}
 463 r'&\rightarrow&r_{\alpha n}=|\vec{r}'-\vec{\tau}_{\alpha n}|\\
 464 \Omega_{\vec{r}'}&\rightarrow&\Omega_{\vec{r'}-\vec{\tau}_{\alpha n}}\\
 465 \delta V_l(r')&\rightarrow&\delta V_l(|\vec{r}'-\vec{\tau}_{\alpha n}|)\\
 466 u_l(r')&\rightarrow&u_l(|\vec{r}'-\vec{\tau}_{\alpha n}|)\\
 467 Y_{lm}(\Omega_{\vec{r}'})&\rightarrow&
 468 Y_{lm}(\Omega_{\vec{r}'-\vec{\tau}_{\alpha n}})
 469 \text{ .}
 470 \end{eqnarray}
 471 Within an iterative treatment on a real space grid consisting of $n_{\text{g}}$ grid points, the sum
 472 \begin{equation}
 473 \sum_{\vec{r}''_{\alpha n}}
 474 \sum_{lm}-i\hbar(\vec{r}'_{\alpha n}\times \nabla_{\vec{r}'_{\alpha n}})
 475 \braket{\vec{r}'_{\alpha n}}{\chi^{\text{SO}}_{lm}}
 476 E^{\text{SO,KB}}_l\braket{\chi^{\text{SO}}_{lm}}{\vec{r}''_{\alpha n}}
 477 \braket{\vec{r}''_{\alpha n}}{\Psi}
 478 \qquad\forall\,\bra{\vec{r}'_{\alpha n}}
 479 \end{equation}
 480 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.
 481 Thus, this method scales linearly with the number of atoms.
 482
 483 The $E_l^{\text{SO,KB}}$ are given by
 484 \begin{equation}
 485 E_l^{\text{SO,KB}}=
 486 \frac{\braket{\delta V_lu_l}{u_l\delta V_l}}
 487      {\bra{u_l}\delta V_l\ket{u_l}}=
 488 \frac{\int_{r}\delta V^2_l(r)u^2_l(r)}r^2dr
 489      {\int_{r'}\int_{r''}\braket{u_l}{r'}\bra{r'}\delta V_l
 490 \ket{r''}\braket{r''}{u_l}}=
 491 \end{equation}
 492