---
abstract: 'We describe the implementation of total angular momentum dependent pseudopotentials in a plane wave formulation of density functional theory. Our approach thus goes beyond the scalar–relativistic approximation usually made in [*ab initio*]{} pseudopotential calculations and explicitly includes spin–orbit coupling. We outline the necessary extensions and compare the results to available all–electron calculations and experimental data.'
author:
- 'G. Theurich'
- 'N. A. Hill'
bibliography:
- '../../PhDThesis/MyBib/mybib.bib'
title: |
    Spin–Orbit Coupling in the\
    [*ab initio*]{} Pseudopotential Framework
---

The [*ab initio*]{} pseudopotential method [@Ihm79p4409; @Payne92p1045; @Kresse96p15] has become a standard tool in many areas of electronic structure calculation. Even magnetic compounds containing $3d$ transition metal ions lie in the realm of the plane wave pseudopotential approach of density functional theory [@Sasaki95p12760; @Moroni97p15629]. In order to obtain high precision results it is necessary to include relativistic effects when calculating the electronic structure of materials containing third row elements [@Bachelet82p2103]. Hence it is now standard procedure to create scalar–relativistic pseudopotentials that include the kinematic relativistic effects (mass–velocity and Darwin term) from the fully relativistic all–electron solution of the atom [@Bachelet82p4199; @Hamann89p2980; @Rappe90p1227; @Troullier91p1993; @Vanderbilt90p7892].

The spin–orbit interaction, however, is only effectively taken into account by the construction of $j$–averaged pseudopotentials for each angular momentum $l$. Thus no spin–orbit splitting is present in the resulting band structure. Although the scalar–relativistic approximation is acceptable in many situations it becomes insufficient in cases where the observed quantities, such as hole effective masses or spin relaxation times, are a direct consequence of the spin–orbit splitting [@Cardona88p1806].

In this paper we report on the implementation of spin–orbit coupling in the pseudopotential scheme. We give the equations required to program the formalism within a generalized spinor approach, and compare the results to experimental data and to all-electron calculations. This is, to our knowledge, the first zeroth order implementation of spin–orbit coupling in the [*ab initio*]{} pseudopotential scheme. Prior publications on this matter always relied on a second variation of the scalar–relativistic zeroth order eigenstates, including spin–orbit coupling to first order in perturbation theory [@Hybertsen86p2920; @Surh91p4286; @Hemstreet93p4238].

Although the fully relativistic treatment of the problem would require a four–current formulation with Dirac spinors it has been shown by Kleinman that a Pauli–like Schrödinger equation captures all relativistic effects to order $\alpha^2$, where $\alpha$ is the fine structure constant [@Kleinman80p2630]. The total ionic pseudopotential to be used is $$\label{j_psp}
V_{PS} =\sum_{l,j,m_j}{{| \Phi^{l,j}_{m_j} \rangle}}V_{l,j}{{\langle \Phi^{l,j}_{m_j} |}}\ ,$$ where the ${{| \Phi^{l,j}_{m_j} \rangle}}$ are the total angular momentum eigenfunctions which can be written in terms of the spherical harmonics, $Y^m_l$, and the eigenfunctions of the $z$–component of the Pauli spin operator, ${{| \uparrow \rangle}}$ and ${{| \downarrow \rangle}}$. For $j=l+\frac{1}{2}$, $m_j=m+\frac{1}{2}$ the ${{| \Phi^{l,j}_{m_j} \rangle}}$ equal $$\left(\frac{l+m+1}{2l+1}\right)^\frac{1}{2}{{| Y^m_l \rangle}}{{| \uparrow \rangle}}+         \left(\frac{l-m}{2l+1}\right)^\frac{1}{2}{{| Y^{m+1}_l \rangle}}{{| \downarrow \rangle}}$$ and for $j=l-\frac{1}{2}$, $m_j=m-\frac{1}{2}$ have the form $$\left(\frac{l-m+1}{2l+1}\right)^\frac{1}{2}{{| Y^{m-1}_l \rangle}}{{| \uparrow \rangle}}-        \left(\frac{l+m}{2l+1}\right)^\frac{1}{2}{{| Y^m_l \rangle}}{{| \downarrow \rangle}}
\label{j_minus}\ .$$ Hence the operator $V_{PS}$ acts in both orbital and spin space. Note that there is only one radial pseudopotential component $V_{l,j}$ with $j=\frac{1}{2}$ for $l=0$ but two with $j=l+\frac{1}{2}$ and $j=l-\frac{1}{2}$ for each $l>0$. The index $m_j$ in equation (\[j\_psp\]) runs from $-j$ to $+j$. It is computationally more efficient to transcribe each term of the semi–local pseudopotential operator $V_{PS}$ into the fully separable Kleinman–Bylander (KB) form [@Kleinman82p1425] $$\label{relKB} V_{KB} = \sum_{i_s,
i_a}\sum_{l,j,m_j}\frac{{{| \delta
V^{i_s,i_a}_{l,j}\phi^{i_s,i_a}_{l,j,m_j} \rangle}}{{\langle \phi^{i_s,i_a}_{l,j,m_j}\delta
V^{i_s,i_a}_{l,j} |}}}{{{\langle \phi^{i_s}_{l,j,m_j} |}}\delta
V^{i_s}_{l,j}{{| \phi^{i_s}_{l,j,m_j} \rangle}}}$$ using the solutions of the atomistic pseudopotential problem $${{| \phi^{i_s,i_a}_{l,j,m_j} \rangle}}={{| R^{i_s,i_a}_{l,j} \rangle}}{{| \Phi^{l,j}_{m_j} \rangle}}\ ,$$ where ${{| R^{i_s,i_a}_{l,j} \rangle}}$ is the radial part of the pseudo eigenfunction of atom species $i_s$ at position $r_{i_s,i_a}$. The potential $\delta V^{i_s,i_a}_{l,j}$ is defined as the difference $$\label{minus_local}
\delta V^{i_s,i_a}_{l,j}(r) = V_{l,j}(r-r_{i_s,i_a}) - V_{loc}(r-r_{i_s,i_a})\ ,$$ where $V_{loc}(r)$ is an arbitrary local potential that needs to be chosen such that the remaining $\delta V$’s are short ranged. The complete KB pseudopotential operator is thus given as the sum of the local part and the non–local KB operator.

To our knowledge all previous pseudopotential calculations that included spin–orbit coupling did so by using a second variation step on the scalar–relativistic zeroth order wave functions, thus including the spin–orbit term to first order perturbation theory [@Hybertsen86p2920; @Surh91p4286; @Hemstreet93p4238]. In contrast we solve directly for general two–component spinor Bloch wave functions expanding in a plane wave spinor basis $$\label{plane_wave}
{{| \psi_{nk} \rangle}} = \sum_{G,\sigma}c^{n,k}_{G,\sigma}{{| k+G \rangle}}{{| \sigma \rangle}},$$ where $G$ are reciprocal lattice vectors and $\sigma$ sums over up and down spin. In the basis of equation (\[plane\_wave\]) the action of the KB operator is as follows:

$$\label{KBelement}
{{\langle \sigma |}}{{\langle k+G |}}V_{KB}{{| \psi_{nk} \rangle}}=\sum_{i_s, i_a}\sum_{l,j,m_j}D^{KB}_{i_s,l,j}\ \varphi^{i_s,i_a}_{k+G}\
M^{KB, \sigma}_{i_s,l,j,m_j,k+G}\ f^{KB}_{i_s,i_a,k,n,l,j,m_j}\ ,$$

where $$D^{KB}_{i_s,l,j} = \left(\frac{4\pi}{V}\right)\frac{1}{\int dr\ r^2 R^{*,i_s}_{l,j}(r)\ \delta V^{i_s}_{l,j}(r)\ R^{i_s}_{l,j}(r)}$$ and $$\varphi^{i_s,i_a}_{k+G} = e^{-i(\vec k+\vec G)\cdot {\vec r}_{i_s,
i_a}}$$ is a phase factor associated with the atomic position. The spin dependent factor $M^{KB}$ of equation (\[KBelement\]) can be written as a spinor and for $j=l+\frac{1}{2}$ is $$\label{M_plus}
M^{KB, \sigma}_{i_s,l,j,m_j,k+G} =
{
\sqrt{l+m+1}\ F^{KB}_{i_s,l,j,m,k+G}
\choose
\sqrt{l-m}\ F^{KB}_{i_s,l,j,m+1,k+G}
}$$ and $$\label{M_minus}
M^{KB, \sigma}_{i_s,l,j,m_j,k+G} =
{
\sqrt{l-m+1}\ F^{KB}_{i_s,l,j,m-1,k+G}
\choose
-\sqrt{l+m}\ F^{KB}_{i_s,l,j,m,k+G}
}$$ for $j=l-\frac{1}{2}$. Also the last factor of equation (\[KBelement\]), $f^{KB}_{i_s,i_a,k,n,l,j,m_j}$, depends on $j$ as follows $$\begin{aligned}
&\sum_{G'} \varphi^{* i_s,i_a}_{k+G}&
(c^{n,k}_{G',\uparrow}\sqrt{l+m+1}\ F^{* KB}_{i_s,l,j,m,k+G'}\nonumber\\
&&+c^{n,k}_{G',\downarrow}\sqrt{l-m}\ F^{* KB}_{i_s,l,j,m+1,k+G'})\end{aligned}$$ for $j=l+\frac{1}{2}$ and $$\begin{aligned}
&\sum_{G'} \varphi^{* i_s,i_a}_{k+G}&
(c^{n,k}_{G',\uparrow}\sqrt{l-m+1}\ F^{* KB}_{i_s,l,j,m-1,k+G'}\nonumber\\
&&-c^{n,k}_{G',\downarrow}\sqrt{l+m}\ F^{* KB}_{i_s,l,j,m+1,k+G'})\end{aligned}$$ for $j=l-\frac{1}{2}$. Finally the KB factors $F^{KB}$ appearing in equations (\[M\_plus\]) and (\[M\_minus\]) are defined as $$\begin{aligned}
F^{KB}_{i_s,l,j,m,k+G} =\sqrt{\frac{4\pi}{2l+1}}Y^m_l(\theta, \varphi)\times\nonumber\\
\int dr\ r^2\ j_l(|k+G|r)\ \delta V^{i_s}_{l,j}(r)\ R^{i_s}_{l,j}(r),\end{aligned}$$ where $Y^m_l$ are the spherical harmonics, the polar angles $\theta$ and $\varphi$ are determined by the vector $\vec k + \vec G$ and $j_l$ are the spherical Bessel functions. The KB factors are calculated once and stored in memory. The contribution of a state ${{| \psi_{nk} \rangle}}$ to the non–local KB part of the total energy is thus given by $$E^{KB}=\sum_{i_s, i_a}\sum_{l,j,m_j}D^{KB}_{i_s,l,j}\left|
f^{KB}_{i_s,i_a,k,n,l,j,m_j}
\right|^2\ .$$

In order to test our implementation we calculated the properties of GaAs and compared with all–electron calculations and available experimental data. Self–consistency was achieved by direct minimization of the total energy via a conjugate gradient method [@Payne92p1045]. The gallium and arsenic pseudopotentials were created following the Troullier–Martins scheme [@Troullier91p1993], and both contained $s$ and $p$ components. Care must be given to the local part of the pseudopotential entering in equation (\[minus\_local\]) to ensure good transferability. We used the $j$–average of the unbound $4d$ state in case of gallium and likewise the $j$–average of the $p$ states for the local part of the arsenic pseudopotential.

                 S-FKKR$^a$   S-FLAPW$^b$   R-PWPP   Exp.$^c$
  ------------- ------------ ------------- -------- ----------
  $a_0$ (Å)         5.56         5.620      5.642     5.653
  $B_0$ (GPa)        77           74         72.2      74.8

  : Equilibrium lattice constant and bulk modulus determined in this work (R-PWPP) compared to all–electron calculations $^a$ Scalar relativistic FKKR, M. Asato et al., PRB [**60**]{}, 5202 (1999), $^b$ Scalar relativistic FLAPW, C. Filippi et al., PRB [**50**]{}, 14947 (1994) and experiment $^c$ Landolt–Börnstein, Vol 22 (1987)[]{data-label="tab_structural"}

The results for lattice constant and bulk modulus for GaAs are shown in table \[tab\_structural\]. The good agreement between our fully relativistic pseudopotential results and the scalar–relativistic all–electron values in table \[tab\_structural\] confirms our approach and reaffirms the notion that spin–orbit splittings have little effect on the structural properties of semiconductors [@Bachelet85p879]. Our calculated bandstructure in figure \[fig\_band\] on the other hand shows clear evidence of spin–orbit coupling. The top of the valence band splits into the light hole, heavy hole manifold and the split off band, separated by 350meV. A similar split is also observed in the upper conduction bands at the Brillouin zone center. In table \[tab\_splitting\] we compare the characteristic spin–orbit splittings we obtained for GaAs at the experimental lattice constant with values from two fully relativistic all–electron calculations found in the literature. The agreement with both all–electron calculations is excellent.

  Splitting                      R-FKKR$^a$   R-FLAPW$^b$   R-PWPP
  ---------------------------- ------------ ------------- --------
  $\Delta_0(\Gamma_{15}^v)$            0.35          0.34     0.35
  $\Delta'_0(\Gamma_{15}^c)$           0.20             –     0.20
  $\Delta_1(L_{3}^v)$                  0.09          0.09     0.09
  $\Delta(X_{5}^v)$                    0.21          0.20     0.22

  : Spin–orbit splittings for GaAs obtained in this work (R-PWPP) compared to the results of two relativistic all–electron calculations: $^a$ Fully relativistic FKKR, S. Bei der Kellen, A. J. Freeman, PRB [**54**]{}, 11187 (1996) and $^b$ Scalar relativistic FLAPW + 2nd variation, C. Filippi, et al., PRB [**50**]{}, 14947 (1994)[]{data-label="tab_splitting"}

For completeness we compare in table \[tab\_eigenvalues\] the eigenvalue spectrum at three special $k$ points of the Brillouin zone with the eigenvalues obtained from the two fully relativistic all–electron calculations cited in table \[tab\_splitting\]. The zero of energy was chosen to coincide with the top of the valence band. Despite the generally good agreement there are two obvious discrepancies at the Brillouin zone center that need some clarification. First the direct band gap of our pseudopotential calculation at $\Gamma$ is more than 4 times larger than the gap resulting from the all–electron calculations. Second the valence band width of our approach is slightly smaller compared to the all–electron results. Both discrepancies result from the fact that the gallium $3d$ orbitals were placed in the frozen core in our calculation but are free to change in the all–electron approaches. Due to the well known self–interaction problem of the local density approximation to density functional theory [@Perdew81p5048] these fairly localized states will lie too high in energy when not frozen. The symmetry of the $d$ states in the zincblende lattice at $\Gamma$ only allows hybridization with $p$ states, e.g. the top of the valence band. Hence the top of the valence band will shift upwards, leading to a reduced band gap and at the same time an increase in the valence band width. Due to the mixed character of the band states away from the Brillouin zone center the effect of the gallium $3d$ states is most pronounced at $\Gamma$.

[c|rr|r]{} Level & R-FKKR$^a$ & R-FLAPW$^b$ & R-PWPP\
$\Gamma_6^v$ & -12.94 & -12.91 & -12.67\
$\Gamma_7^v$ & -0.35 & -0.34 & -0.35\
$\Gamma_8^v$ & 0.00 & 0.00 & 0.00\
$\Gamma_6^c$ & 0.12 & 0.17 & 0.69\
$\Gamma_7^c$ & 3.46 & & 3.40\
$\Gamma_8^c$ & 3.66 & & 3.60\
\
$X_6^v$ & -10.42 & -10.41 & -10.35\
$X_7^v$ & -7.02 & -7.00 & -6.83\
$X_6^v$ & -2.88 & -2.85 & -2.74\
$X_7^v$ & -2.79 & -2.76 & -2.65\
$X_6^c$ & 1.17 & 1.23 & 1.23\
$X_7^c$ & 1.39 & & 1.44\
\
$L_6^v$ & -11.18 & -11.14 & -11.06\
$L_6^v$ & -6.83 & -6.82 & -6.63\
$L_6^v$ & -1.38 & -1.36 & -1.32\
$L_{4,5}^v$ & -1.17 & -1.16 & -1.10\
$L_6^c$ & 0.71 & 0.73 & 0.97\
$L_6^c$ & 4.38 & & 4.34\
$L_{4,5}^c$ & 4.46 & & 4.44\

Compared to calculations that do not include the spin–orbit term we find that the inclusion of spin–orbit coupling worsens the short coming of the local density approximation of underestimating the band gap. The reason for this observation simply lies in the fact that the top of the valence band splits and the light and heavy hole states move closer to the bottom of the conduction band.

In conclusion, we have implemented spin–orbit coupling in the well established [*ab initio*]{} pseudopotential approach of density functional theory. This paper gives the necessary expressions in a two–component spinor plane wave basis and demonstrates the applicability of the method for bulk GaAs. Our results compare very well to relativistic all–electron calculations. Since our direct approach is based on a complete spinor plane wave basis it can easily be extended to systems that show exchange splitting and exhibit non–collinear spin arrangements.

The code will be available under the GNU General Public License [@GPL] at <http://www.mrl.ucsb.edu/~theurich/Spinor/>.

This work was supported by the ONR grant number N00014-00-10557, by NSF-DMR under the grant 9973076 and by ACS PRF under the grant 33851-G5. G.T. acknowledges fellowship support from the UCSB Materials Research Lab., funded by the Corning Foundation.
