text:
# Consequences of local gauge symmetry in empirical tight-binding theory

## I Introduction

Tight-binding theory was originally proposed as an ab initio technique for calculating the electronic properties of crystalline solids from atomic wave functions.Bloch (1928) However, first-principles calculations based on a linear combination of atomic orbitals (LCAO) are computationally very demanding, and the tight-binding approach met with relatively little success until Slater and Koster suggested that it be used as an interpolation scheme,Slater and Koster (1954) in which the Hamiltonian matrix elements are fitted to experimental data or to band structures computed by other methods. This made it possible to describe atomic-level physics in a basis of minimal size, leading to wide-ranging applications in many areas of condensed-matter physics. Harrison (1989, 1999); Wang and Ho (1996); Bullett (1980); Goringe et al. (1997); Turchi et al. (1998) With modern computer capabilities, first-principles electronic-structure calculations are now commonplace, and ab initio tight-binding theories are flourishing. Bullett (1980); Goringe et al. (1997); Turchi et al. (1998) Yet even today, the empirical theory Slater and Koster (1954) predominates (even for the fitting of first-principles calculations) because it is simple and physically intuitive.

The formalism of Slater and Koster (1954) is incomplete, however, in that it contains no prescription for coupling the electronic system to external electromagnetic fields. In ab initio theories, Bullett (1980); Goringe et al. (1997); Turchi et al. (1998) one can use minimal coupling (with suitable modifications for nonlocal potentials Ismail-Beigi et al. (2001)) and calculate directly the necessary matrix elements of the momentum or velocity operator. In the empirical theory, these matrix elements can simply be treated as extra fitting parameters, Chang and Schulman (1985); Chang and Aspnes (1990); Xu (1990) determined by fitting the dielectric function (and thus oscillator strengths) to experimental or first-principles spectra. However, even with the full use of symmetry restrictions, the number of additional parameters can be undesirably large; for example, Chang and Aspnes (1990) have proposed an $`sp^3d^2`$ model for GaAs with 13 Hamiltonian parameters and 17 independent momentum parameters.

It is therefore clearly desirable to find ways of reducing or eliminating these extra fitting parameters. One possibility is to define a kinematic momentum operator (equal to mass $`m`$ times velocity) by

$$𝐩=\frac{m}{i\mathrm{}}[𝐱,H],$$
(1)
where x is the coordinate of the electron and $`H`$ is the Hamiltonian. In a sense this merely trades one problem for another, since the coordinate matrix elements are still unknown, and the number of these allowed by symmetry is no less than the number of momentum matrix elements. However, it is physically reasonable to simplify the coordinate matrix by setting all nonlocal matrix elements to zero:

$$\alpha ,𝐱_i|𝐱|\alpha ^{},𝐱_i^{}=\delta _{ii^{}}[\delta _{\alpha \alpha ^{}}𝐱_i+𝐱_{\alpha \alpha ^{}}(i)].$$
(2)
Here $`|\alpha ,𝐱_i`$ is the ket vector for an orthogonalized atomic orbital (Löwdin orbital Löwdin (1950); Altmann (1991)) of type $`\alpha `$ located at position $`𝐱_i`$. The parameter $`𝐱_{\alpha \alpha ^{}}`$ is an intra-atomic matrix element coupling orbitals $`\alpha `$ and $`\alpha ^{}`$ on the same atom.

The simplest choice of all is to set

$$𝐱_{\alpha \alpha ^{}}0;$$
(3)
in this model, there are no fitting parameters beyond those found in the Hamiltonian. Dresselhaus and Dresselhaus (1967); Smith (1979); Brey and Tejedor (1983); Harrison (1989, 1999); Koiller et al. (1991); Tserbak et al. (1993); Lew Yan Voon and Ram-Mohan (1993); Boykin (1999) A closely related approach is the Peierls substitution, Peierls (1933); Graf and Vogl (1995); Dumitrică et al. (1998); Boykin et al. (2001); Ismail-Beigi et al. (2001); Boykin and Vogl (2002) in which the zero-field Hamiltonian matrix $`\alpha ,𝐱_i|H|\alpha ^{},𝐱_i^{}`$ for a particle of charge $`e`$ is replaced by

$$\alpha ,𝐱_i|H|\alpha ^{},𝐱_i^{}\mathrm{exp}\left(\frac{ie}{\mathrm{}c}_{𝐱_i^{}}^{𝐱_i}𝐀\text{d}𝐱\right)+e\varphi (𝐱_i)\delta _{ii^{}}\delta _{\alpha \alpha ^{}}$$
(4)
in the presence of a vector potential $`𝐀`$ and scalar potential $`\varphi `$. If the path of integration in Eq. (4) is chosen to be a straight line, Graf and Vogl (1995); Boykin et al. (2001); Ismail-Beigi et al. (2001) then the linear term in the Taylor series expansion of this equation is the same as the $`𝐀𝐩`$ coupling obtained from Eqs. (1)–(3). Graf and Vogl (1995)

The total elimination of extra fitting parameters makes this model an attractive one. However, by eliminating the intra-atomic matrix elements $`𝐱_{\alpha \alpha ^{}}`$, one obtains a tight-binding model that is not valid in the tight-binding limit of isolated atoms. Thus, although the model should provide a reasonable description of inter-atomic transitions between extended states, one has less confidence in its ability to describe localized states, which may be important at surfaces or interfaces. Many authors have therefore suggested augmenting the zero-parameter model by including a small number of intra-atomic matrix elements. Selloni et al. (1986); Reining et al. (1994); Bennetto and Vanderbilt (1996); Leung and Whaley (1997); Leung et al. (1998); Cruz et al. (1999); Pedersen et al. (2001) It has been shown for porous Si that, although the intra-atomic matrix elements are small in magnitude (in a Bloch-function basis Bloch (1928)), the interference between these terms and the inter-atomic matrix elements contributes 25% of the total absorption. Cruz et al. (1999) Thus, it appears that a quantitative treatment of nanostructures may not be possible (in general) without the inclusion of intra-atomic matrix elements.

The main difficulty with such models Selloni et al. (1986); Reining et al. (1994); Bennetto and Vanderbilt (1996); Leung and Whaley (1997); Leung et al. (1998); Cruz et al. (1999); Pedersen et al. (2001) is that they are not gauge invariant. Boykin and Vogl (2002) As shown by the examples in Refs. Ismail-Beigi et al., 2001 and Governale and Ungarelli, 1998, lack of gauge invariance can lead to gross qualitative errors in the predicted values of physical quantities. Thus, there are significant problems with both approaches considered above. The models with $`𝐱_{\alpha \alpha ^{}}=0`$ are gauge invariant, but they cannot describe intra-atomic transitions. The models with $`𝐱_{\alpha \alpha ^{}}0`$ can describe intra-atomic transitions, but they are not gauge invariant.

The purpose of this paper is to demonstrate a technique for constructing tight-binding models that are gauge invariant and provide a full description of intra-atomic transitions. This is achieved by treating empirical tight-binding theory not as an approximation derived from the Schrödinger equation, but as a fundamental quantum-mechanical system in its own right. This theory is required to satisfy all of the basic principles of quantum mechanics, the most important of which (in the present context) is the principle of local gauge symmetry. Weyl (1929); Dirac (1931); Yang and Mills (1954); Wu and Yang (1975); Moriyasu (1983); Guidry (1991) The essence of this principle is the concept that electromagnetism in quantum mechanics is the gauge-invariant manifestation of a nonintegrable (i.e., path-dependent) phase factor. Dirac (1931); Wu and Yang (1975)

As will be shown in Sec. II below, the reason why existing models with intra-atomic coupling Selloni et al. (1986); Reining et al. (1994); Bennetto and Vanderbilt (1996); Leung and Whaley (1997); Leung et al. (1998); Cruz et al. (1999); Pedersen et al. (2001) are not gauge invariant is that the coordinates $`x`$, $`y`$, and $`z`$ do not commute. Boykin and Vogl (2002) Gauge invariance requires commuting coordinates, the existence of which implies the existence of a basis of coordinate eigenkets. Since empirical tight-binding theory deals with finite vector spaces, the coordinate basis is necessarily discrete. Hence, the most general gauge-invariant finite vector space is a set of discrete coordinate eigenkets. This basis may be transformed to a tight-binding basis by constructing “orbitals” from symmetrized combinations of coordinate eigenkets (using well-known techniques for symmetrizing plane waves Tinkham (1964); Bassani (1966)).

The concept of gauge symmetry on a discrete lattice is not new, having appeared many years ago as a technique for imposing a momentum cutoff in quantum chromodynamics. Guidry (1991); Wilson (1974); Kogut and Susskind (1975); Rebbi (1983); Rothe (1997) Governale and Ungarelli (1998) have recently suggested using lattice gauge techniques in empirical tight-binding theory. However, their proposal, like most applications of lattice gauge theory, is based on a simple cubic lattice. As shown below, the simple cubic lattice is unsuitable for practical tight-binding models because it can only achieve sufficient accuracy with an unreasonably large basis (i.e., a very small lattice constant). Thus, the development of efficient tight-binding models requires consideration of more general geometries.

Christ, Friedberg, and Lee (1982a, b, c) have developed a lattice gauge theory for random lattices, which (with some slight modifications) is sufficiently general for the present purposes. However, the complete formal machinery of quantum chromodynamics is somewhat cumbersome when one is dealing only with simple electromagnetism. Thus, for reasons of clarity, the author has chosen to present the theory in terms of a simple but elegant approach used by Dirac. Dirac (1931) After a preliminary discussion of topology (i.e., how an electron is permitted to move from one lattice site to another) in Sec. III, Sec. IV presents an adaptation of Dirac’s analysis Dirac (1931) to the case of a discrete lattice. The outcome is a gauge-invariant formulation of electromagnetism in empirical tight-binding theory.

Although the theory derived in this way has many similarities with conventional tight-binding theory, there are significant differences as well. Not all tight-binding models can be made gauge invariant; not (a) this is possible only if the basis can be constructed from symmetrized coordinate eigenkets. In addition, local gauge symmetry imposes strong restrictions on the Hamiltonian matrix, which have the effect of sharply reducing the number of allowed Hamiltonian fitting parameters. Finally, unlike previous empirical tight-binding theories, the present approach provides an explicit (discrete) wave function for the electron.

The formalism derived here is applied to two semiconductors with the diamond structure (Ge and Si) in Sec. V. For these systems, a basis of 15 orbitals per atom is shown to provide a satisfactory fit to the valence bands and the lowest conduction bands (up to about 5 eV above the valence-band maximum). These results are comparable to those obtained from a 10-orbital basis proposed recently in Ref. Jancu et al., 1998. The basis used here is 50% larger, but their model Jancu et al. (1998) cannot be made gauge invariant if intra-atomic coupling is included. Thus, it appears that some tradeoffs are necessary if gauge invariance is to be achieved.

## II Coordinate matrices and the coordinate representation

As mentioned above, the intra-atomic coupling used in existing tight-binding models Selloni et al. (1986); Reining et al. (1994); Bennetto and Vanderbilt (1996); Leung and Whaley (1997); Leung et al. (1998); Cruz et al. (1999); Pedersen et al. (2001) leads to a lack of gauge invariance. Boykin and Vogl (2002) This may be seen from a simple $`sp^3`$ model for a single atom. In this case, we know from atomic physics that there are coordinate matrix elements coupling the $`s`$ and $`p`$ orbitals:

$$s|x|p_x=s|y|p_y=s|z|p_zc,$$
(5)
where $`c`$ is real. In the basis $`\{|s,|p_x,|p_y,|p_z\}`$, the matrices representing $`x`$ and $`y`$ are therefore

$$x=\left[\begin{array}{cccc}0& c& 0& 0\\ c& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],y=\left[\begin{array}{cccc}0& 0& c& 0\\ 0& 0& 0& 0\\ c& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].$$
(6)
But this implies that $`x`$ and $`y`$ do not commute:

$$xyyx=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& c^2& 0\\ 0& c^2& 0& 0\\ 0& 0& 0& 0\end{array}\right]0.$$
(7)
This means that the coordinate representation (consisting of simultaneous eigenkets of $`x`$, $`y`$, and $`z`$) does not exist. Even more important, it means that the theory cannot be gauge invariant. In a gauge transformation, the vector and scalar potentials transform as

$$𝐀𝐀+\mathrm{\Lambda },\varphi \varphi \frac{1}{c}\frac{\mathrm{\Lambda }}{t},$$
(8)
and the state ket $`|\psi `$ transforms as

$$|\psi U|\psi ,U\mathrm{exp}\left[\frac{ie\mathrm{\Lambda }(𝐱,t)}{\mathrm{}c}\right],$$
(9)
where $`\mathrm{\Lambda }`$ is an arbitrary function of $`𝐱=(x,y,z)`$ and $`t`$. If a theory is gauge invariant, all physically measurable quantities must be independent of such transformations. But the expectation value $`x`$ is a measurable quantity, and under a gauge transformation one has

$$xU^{}xU,$$
(10)
where $`U^{}xUx`$ if $`\mathrm{\Lambda }`$ depends on $`y`$ or $`z`$. Hence, no theory can be gauge invariant if $`x`$, $`y`$, and $`z`$ do not commute.

Is there any way of achieving gauge invariance without setting $`c=0`$? Perhaps the $`sp^3`$ basis is too small, and the situation might be improved by including more orbitals ($`d,f,\mathrm{}`$). However, one soon finds that for any finite LCAO basis, the lack of gauge invariance persists. This follows directly from the Wigner-Eckart theorem—since $`𝐱`$ is a vector operator, it couples states with angular momentum $`l`$ to those with $`l\pm 1`$. Hence, any finite truncation of the basis results in non-commuting coordinates.

Another possibility is to keep the same $`sp^3`$ basis, but modify the coordinate matrix. The physical justification for doing so is the fact that the orbitals used in empirical tight-binding theory are not atomic orbitals; they are orthogonalized atomic orbitals. Löwdin (1950) Therefore, they do not have the full rotational symmetry of atomic orbitals; they have only the site symmetry of the crystal structure. For example, the atoms in a diamond crystal have site symmetry $`T_d`$. Hahn (1992) Therefore, the orbital that was denoted $`|p_z`$ above should really be written as $`|\mathrm{\Gamma }_{15}^z`$, since it belongs to the $`\mathrm{\Gamma }_{15}`$ representation of $`T_d`$. Parmenter (1955)

However, the $`d`$ orbital $`|d_{xy}`$ also transforms as $`|\mathrm{\Gamma }_{15}^z`$. Thus, in the $`T_d`$ group, the matrix element $`b\mathrm{\Gamma }_{15}^x|y|\mathrm{\Gamma }_{15}^z`$ is allowed, and the coordinate matrices (6) become

$$x=\left[\begin{array}{cccc}0& c& 0& 0\\ c& 0& 0& 0\\ 0& 0& 0& b\\ 0& 0& b& 0\end{array}\right],y=\left[\begin{array}{cccc}0& 0& c& 0\\ 0& 0& 0& b\\ c& 0& 0& 0\\ 0& b& 0& 0\end{array}\right].$$
(11)
This yields

$$xyyx=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& c^2b^2& 0\\ 0& b^2c^2& 0& 0\\ 0& 0& 0& 0\end{array}\right],$$
(12)
which is equal to zero if $`b=\pm c`$.

Setting $`b=c`$ would imply that the $`p`$ and $`d`$ orbitals have equal weight in the $`\mathrm{\Gamma }_{15}`$ states. This is not as absurd as it sounds; Boguslawski and Gorczyca (1994) have shown using first-principles pseudopotential calculations that for the $`\mathrm{\Gamma }_{15}`$ states at the top of the valence band in GaAs, the probability of finding an electron in a cation $`d`$ orbital is greater than that of finding it in a cation $`p`$ orbital. (In AlAs, the probability ratio is greater than 2.Boguslawski and Gorczyca (1994)) Thus, it is not unreasonable to assume that $`b`$ and $`c`$ have comparable magnitudes.

If one sets $`b=c`$, then the coordinate operators have simultaneous eigenkets $`|x^{},y^{},z^{}`$, which are given by

$`|ccc`$ $`=`$ $`\frac{1}{2}(|\mathrm{\Gamma }_1+|\mathrm{\Gamma }_{15}^x+|\mathrm{\Gamma }_{15}^y+|\mathrm{\Gamma }_{15}^z),`$
$`|c\overline{c}\overline{c}`$ $`=`$ $`\frac{1}{2}(|\mathrm{\Gamma }_1+|\mathrm{\Gamma }_{15}^x|\mathrm{\Gamma }_{15}^y|\mathrm{\Gamma }_{15}^z),`$
$`|\overline{c}c\overline{c}`$ $`=`$ $`\frac{1}{2}(|\mathrm{\Gamma }_1|\mathrm{\Gamma }_{15}^x+|\mathrm{\Gamma }_{15}^y|\mathrm{\Gamma }_{15}^z),`$
$`|\overline{c}\overline{c}c`$ $`=`$ $`\frac{1}{2}(|\mathrm{\Gamma }_1|\mathrm{\Gamma }_{15}^x|\mathrm{\Gamma }_{15}^y+|\mathrm{\Gamma }_{15}^z).`$ (13)
Note that the coordinate eigenvalues are located at the corners of a tetrahedron. In fact, the linear combinations given in (13) are identical to the hybrid bond orbitals used in analytical tight-binding theories, Harrison (1989, 1999) although these are not ordinarily interpreted as exact coordinate eigenkets because the $`\mathrm{\Gamma }_{15}`$ states are assumed to be pure $`p`$ orbitals.

The procedure outlined above is rather clumsy; one simply modifies the coordinate matrices by trial and error in an attempt to make them commute. One cannot predict in advance whether the attempt will succeed, and in general it will not. However, the unitary transformation (13) may be inverted to obtain

$`|\mathrm{\Gamma }_1`$ $`=`$ $`\frac{1}{2}(|ccc+|c\overline{c}\overline{c}+|\overline{c}c\overline{c}+|\overline{c}\overline{c}c),`$
$`|\mathrm{\Gamma }_{15}^x`$ $`=`$ $`\frac{1}{2}(|ccc+|c\overline{c}\overline{c}|\overline{c}c\overline{c}|\overline{c}\overline{c}c),`$
$`|\mathrm{\Gamma }_{15}^y`$ $`=`$ $`\frac{1}{2}(|ccc|c\overline{c}\overline{c}+|\overline{c}c\overline{c}|\overline{c}\overline{c}c),`$
$`|\mathrm{\Gamma }_{15}^z`$ $`=`$ $`\frac{1}{2}(|ccc|c\overline{c}\overline{c}|\overline{c}c\overline{c}+|\overline{c}\overline{c}c),`$ (14)
which immediately suggests a more fruitful approach. The linear combinations given in Eq. (14) are just what one would obtain by starting with a single coordinate eigenket (say $`|ccc`$) and using the symmetry operations of the tetrahedral group $`T_d`$ to construct “symmetrized” orbitals Tinkham (1964); Bassani (1966) that transform according to the irreducible representations of $`T_d`$. Parmenter (1955)

Thus, in this alternative approach, the coordinate basis is taken as fundamental, and the tight-binding basis is merely a secondary alternative that is useful for reasons of symmetry. Since the existence of a coordinate representation is necessary for gauge invariance, no tight-binding basis can be made gauge invariant if it cannot be represented in terms of symmetrized coordinate eigenkets. not (a) Hence, this symmetrization procedure provides us with all possible gauge-invariant tight-binding models.

The orbitals in Eq. (14) are useful as a starting point, but they cannot be interpreted as atomic orbitals because they do not have inversion symmetry. To obtain more “atomic-like” orbitals, one can apply the symmetry operations of the cubic group $`O_h`$ to the basis ket $`|ccc`$, which yields the orbitals

$`|\mathrm{\Gamma }_1`$ $`=`$ $`|s`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{8}}}(|ccc+|c\overline{c}\overline{c}+|\overline{c}c\overline{c}+|\overline{c}\overline{c}c`$
$`+|\overline{c}\overline{c}\overline{c}+|\overline{c}cc+|c\overline{c}c+|cc\overline{c}),`$
$`|\mathrm{\Gamma }_2^{}`$ $`=`$ $`|f_{xyz}`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{8}}}(|ccc+|c\overline{c}\overline{c}+|\overline{c}c\overline{c}+|\overline{c}\overline{c}c`$
$`|\overline{c}\overline{c}\overline{c}|\overline{c}cc|c\overline{c}c|cc\overline{c}),`$
$`|\mathrm{\Gamma }_{15}^z`$ $`=`$ $`|p_z`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{8}}}(|ccc|c\overline{c}\overline{c}|\overline{c}c\overline{c}+|\overline{c}\overline{c}c`$
$`|\overline{c}\overline{c}\overline{c}+|\overline{c}cc+|c\overline{c}c|cc\overline{c}),`$
$`|\mathrm{\Gamma }_{25^{}}^{xy}`$ $`=`$ $`|d_{xy}`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{8}}}(|ccc|c\overline{c}\overline{c}|\overline{c}c\overline{c}+|\overline{c}\overline{c}c`$
$`+|\overline{c}\overline{c}\overline{c}|\overline{c}cc|c\overline{c}c+|cc\overline{c}).`$
Here two different labels are used: the representations of $`O_h`$, Bouckaert et al. (1936) and the conventional atomic orbital notation. Other orbitals not given here may be obtained from cyclic permutations of $`x`$, $`y`$, and $`z`$. Note that the $`T_d`$ orbital $`|\mathrm{\Gamma }_{15}^z`$ in Eq. (14) is the same as $`\frac{1}{\sqrt{2}}(|p_z+|d_{xy})`$, whereas the $`T_d`$ orbital $`|\mathrm{\Gamma }_1`$ is just $`\frac{1}{\sqrt{2}}(|s+|f_{xyz})`$.

In the basis (II), the nonzero coordinate matrix elements are

$$s|x|p_x=p_x|y|d_{xy}=d_{xy}|z|f_{xyz}=c,$$
(16)
plus others given by cyclic permutations. The selection rules for $`𝐱`$ are thus the same as those in a spherically symmetric atom—although, in any real atom, the matrix elements (16) would not be numerically equal. This equality occurs because the basis kets (II) are degenerate eigenkets of the radial coordinate $`r=\sqrt{x^2+y^2+z^2}`$. To break the numerical equality, one would need to use basis functions with a linear combination of different radii.

The above procedure may, of course, be applied to coordinate eigenkets other than $`|ccc`$. Below is a list of the representations obtained by applying the symmetry operations of $`O_h`$ to several different “generator” eigenkets:

$`|000`$ $``$ $`\mathrm{\Gamma }_1`$
$``$ $`s,`$
$`|100`$ $``$ $`\mathrm{\Gamma }_1+\mathrm{\Gamma }_{15}+\mathrm{\Gamma }_{12}`$
$``$ $`s^1p^3d^2,`$
$`|111`$ $``$ $`\mathrm{\Gamma }_1+\mathrm{\Gamma }_{15}+\mathrm{\Gamma }_{25^{}}+\mathrm{\Gamma }_2^{}`$
$``$ $`s^1p^3d^3f^1,`$
$`|110`$ $``$ $`\mathrm{\Gamma }_1+\mathrm{\Gamma }_{12}+\mathrm{\Gamma }_{15}+\mathrm{\Gamma }_{25^{}}+\mathrm{\Gamma }_{25}`$ (17)
$``$ $`s^1p^3d^5f^3.`$
Explicit basis functions for these representations are given in Appendix A.

With these results, one can now construct a gauge-invariant tight-binding model simply by putting a set (or more than one set) of these “orbitals” on each atom in a crystal or molecule. In such a model, the coordinate operators $`x`$, $`y`$, and $`z`$ commute by construction. However, one is no longer permitted to choose orbitals arbitrarily. The choices are limited to taking all of the orbitals in a set or taking none of them. As an example, one cannot discard the $`f`$ orbitals in the basis generated by $`|110`$ without destroying the gauge invariance of the theory.

This approach yields a tight-binding model with orthogonal orbitals. Another approach is to define a grid of coordinates, some points of which are not uniquely associated with individual atoms. One may still construct symmetrized orbitals in this case, but the orbitals are not orthogonal. This makes the tight-binding approach more difficult; however, one can simplify the theory by choosing a Bravais lattice for the coordinate grid, in which case the model may be viewed as a discrete pseudopotential model. Applications of both the pseudopotential and tight-binding approaches are considered below in Sec. V.

## III Topology of the lattice

As we have seen, the most general gauge-invariant tight-binding basis consists of a set of discrete coordinate eigenkets, which will be referred to as a lattice. Such a lattice is generally not periodic. In order to apply the principle of local gauge symmetry to such a system, one must be able to calculate the change in phase that occurs along any specified path in coordinate space. Dirac (1931); Wu and Yang (1975) Thus, the first step is to define what is meant by a “path” in a discrete coordinate system.

In general, a path is just an ordered sequence of points. In a continuous coordinate system, neighboring points in the sequence must be separated by an infinitesimal distance. This defines the topology of the system, in which points are linked together only if they are adjacent in coordinate space. It is desirable to define the topology of the discrete lattice in a similar way.

One way to do this is to construct a Voronoi polyhedron around each site in the lattice. Christ et al. (1982a) A Voronoi polyhedron is just the region in space closest to that point; Okabe et al. (2000) if the lattice is a Bravais lattice, the polyhedron is the same as a Wigner-Seitz cell. In mathematical terms, the Voronoi polyhedron $`\mathrm{\Omega }_i`$ for site $`𝐱_i`$ is the set of points $`𝐱`$ such that $`|𝐱𝐱_i||𝐱𝐱_j|`$ for all $`ji`$. The topology is then defined by the following rule: If $`\mathrm{\Omega }_i`$ and $`\mathrm{\Omega }_j`$ share a surface with area $`S_{ij}>0`$, the sites $`𝐱_i`$ and $`𝐱_j`$ are linked together; otherwise, they are not.

In some cases, it may happen that two Voronoi polyhedra share only a point or a line, in which case $`S_{ij}=0`$. The linking algorithm presented in Ref. Christ et al., 1982a does not consider this possibility (because Ref. Christ et al., 1982a deals only with random lattices, for which the probability of such an event is zero). For certain applications, it is useful to include links between such sites, Okabe et al. (2000) but we shall see below that these links should be excluded in the present situation. Thus, only adjacent sites whose Voronoi polyhedra share a surface with $`S_{ij}>0`$ are linked.

A path in the discrete lattice is then just an ordered sequence of linked points, and a closed path is one whose first and last points are the same. By definition, every edge of a Voronoi polyhedron is equidistant from three or more lattice sites, all of which lie in a plane perpendicular to the given edge. These sites are closer to this edge than any other sites. The links between these sites form a closed path, and the area bordered by the links is called a “plaquette.” There is a one-to-one relationship between the plaquettes and the edges of the Voronoi polyhedra.

The plaquettes partition all of coordinate space into nonoverlapping volumes called (Delaunay) cells. Each cell is uniquely associated with one corner of a Voronoi polyhedron. The partitioning of space into cells is referred to as a Delaunay tessellation. Okabe et al. (2000)

A general algorithm for calculating the geometry of Voronoi polyhedra, links, plaquettes, and cells is presented in Appendix B. The expressions derived there will be of use in what follows.

## IV Local gauge symmetry on an arbitrary discrete lattice

Christ, Friedberg, and Lee (1982a, b, c) have developed a theory of local gauge symmetry on a random lattice. This section presents a modified version of their theory, with special emphasis on the implications of the principle of local gauge symmetry for tight-binding theory. The presentation follows Dirac’s approach, Dirac (1931) in which the existence of electromagnetic fields is “derived” as a straightforward consequence of a degree of freedom (nonintegrable phases) possessed by any quantum-mechanical system that can be represented in a coordinate basis.

### IV.1 Electromagnetism is a nonintegrable phase

In a discrete coordinate basis, any ket vector may be expressed as

$$|\psi =\underset{i}{}c_i|𝐱_i,$$
(18)
where $`|𝐱_i`$ is a coordinate eigenket, which is assumed to be normalized such that

$$𝐱_i|𝐱_j=\delta _{ij}.$$
(19)
Dirac’s starting point Dirac (1931) is the fact that physical predictions in quantum mechanics are ultimately expressed in terms of probabilities of the form $`|\psi |\psi ^{}|^2`$, where the probability amplitude $`\psi |\psi ^{}`$ is given by

$$\psi |\psi ^{}=\underset{i}{}c_i^{}c_i^{}.$$
(20)
The probability is obviously well defined even when the overall phase of $`|\psi `$ has no definite value. This degree of freedom is referred to as global gauge symmetry.

The existence of global gauge symmetry raises the question of whether it is necessary for the local probability amplitude $`c_i`$ to have a definite phase. In other words, suppose we write

$$c_i=b_ie^{i\beta _i},$$
(21)
where the phase of $`b_i`$ is well defined (to within an integer multiple of $`2\pi `$), but $`\beta _i`$ is a nonintegrable function—that is, the change in $`\beta _i`$ around a closed path can take on any value. In this case one can see that $`|\psi |\psi ^{}|^2`$ is well defined only if the change in $`\beta _i`$ around any closed path is the same for all states $`|\psi `$ and $`|\psi ^{}`$ (to within an integer multiple of $`2\pi `$, which is absorbed into the definition of $`b_i`$). But anything that is the same for all states can be viewed as a physically real part of the dynamical system. Since the present system consists only of a single point particle, these nonintegrable phases must represent a field of force acting on the particle.

The principle of local gauge symmetry is therefore defined by the following two postulates: Dirac (1931) (i) The physical predictions of the theory must be unambiguous. (ii) The phase of $`c_i`$ at any point in space and time need not be well defined; only the change in phase between linked points must be definite. As shown above, these postulates entail that the change in $`\beta _i`$ around any closed path must be the same for all states. According to postulate (ii), this change is fixed for any path by the change in $`\beta _i`$ between two linked points in space:

$$\kappa _{ij}=\beta _i\beta _j(i\text{ linked to }j),$$
(22)
and in time:

$$\lambda _i=\frac{\text{d}\beta _i}{\text{d}t}\dot{\beta }_i.$$
(23)
Since $`\beta _i`$ is nonintegrable, $`\kappa _{ij}`$ and $`\lambda _i`$ are independent variables. These two quantities are the fundamental dynamical variables that arise from the principle of local gauge symmetry. It will now be shown that $`\kappa _{ij}`$ and $`\lambda _i`$ can be interpreted as potentials for the electromagnetic field.

One possible closed path involves a space displacement $`𝐝_{ij}=𝐱_i𝐱_j`$ and an infinitesimal time displacement $`\text{d}t`$, followed by $`𝐝_{ji}`$ and $`\text{d}t`$. The change in phase around this path can be used to define (tentatively) an electric field variable

$$E_{ji}=\frac{\mathrm{}}{e}\frac{\lambda _j\lambda _i\dot{\kappa }_{ji}}{d_{ji}},$$
(24a)
where $`d_{ji}=|𝐝_{ji}|`$. If the index $`\mathrm{}`$ is used to label the links, one may write this in the simpler form
$$E_{\mathrm{}}=\frac{\mathrm{}}{e}\frac{\mathrm{\Delta }\lambda _{\mathrm{}}\dot{\kappa }_{\mathrm{}}}{d_{\mathrm{}}}.$$
(24b)
Here $`E_{\mathrm{}}`$ is interpreted as the component of the electric field in the direction $`𝐝_{\mathrm{}}=𝐝_{ji}`$; the components perpendicular to $`𝐝_{\mathrm{}}`$ are not defined. Equation (24) takes a familiar form if expressed in terms of the potentials

$$\varphi _i=\frac{\mathrm{}}{e}\lambda _i,A_{\mathrm{}}=\frac{\mathrm{}c}{e}\frac{\kappa _{\mathrm{}}}{d_{\mathrm{}}},$$
(25)
since then

$$E_{\mathrm{}}=\frac{\mathrm{\Delta }\varphi _{\mathrm{}}}{d_{\mathrm{}}}\frac{1}{c}\frac{A_{\mathrm{}}}{t}.$$
(26)
Here the notation $`A_{\mathrm{}}/t`$ indicates that $`d_{\mathrm{}}`$ is to be held constant during the differentiation.

Another type of closed path is an elementary plaquette $`q`$ constructed from links $`\mathrm{}`$ in coordinate space (see Sec. III and Appendix B). The change in phase around the perimeter of the plaquette may be used to define the magnetic field

$$B_q=\frac{1}{S_q}\frac{\mathrm{}c}{e}\underset{\mathrm{}q}{}\kappa _{\mathrm{}}=\frac{1}{S_q}\underset{\mathrm{}q}{}A_{\mathrm{}}d_{\mathrm{}},$$
(27)
where $`S_q`$ is the area of the plaquette \[see Eq. (138)\]. Summing Eq. (27) over the (closed) surface of a cell $`c`$ leads immediately to the “no monopoles” law: not (b)

$$\underset{qc}{}B_qS_q=0.$$
(28)
Likewise, summing $`E_{\mathrm{}}d_{\mathrm{}}`$ around the perimeter of a plaquette gives Faraday’s law:

$$\underset{d_{\mathrm{}}S_q}{}E_{\mathrm{}}d_{\mathrm{}}=\frac{1}{c}\frac{\text{d}}{\text{d}t}(B_qS_q).$$
(29)
The other two Maxwell equations can be obtained from the Lagrangian $`L=L_f+L_e`$, where $`L_f`$ is the electromagnetic field Lagrangian

$$L_f=\frac{1}{8\pi }\underset{\mathrm{}}{}3E_{\mathrm{}}^2\mathrm{\Omega }_{\mathrm{}}\frac{1}{8\pi }\underset{q}{}3B_q^2\mathrm{\Omega }_q.$$
(30)
Here $`\mathrm{\Omega }_{\mathrm{}}=\frac{1}{3}S_{\mathrm{}}d_{\mathrm{}}`$ is the volume of link $`\mathrm{}`$, where $`S_{\mathrm{}}=S_{ij}`$ is the area of the surface shared by the Voronoi polyhedra for sites $`i`$ and $`j`$ (see Appendix B). Likewise, $`\mathrm{\Omega }_q=\frac{1}{3}S_qd_q`$ is the volume of plaquette $`q`$, where $`d_q`$ is the length of the Voronoi polyhedron edge corresponding to $`q`$. Equation (30) is just a discrete version of the standard field Lagrangian Jackson (1975) $`\frac{1}{8\pi }(E^2B^2)\text{d}^3x`$; the only apparent difference is an extra factor of 3. This factor cancels the factor of $`\frac{1}{3}`$ in the definition of $`\mathrm{\Omega }_{\mathrm{}}`$ and $`\mathrm{\Omega }_q`$, thus leading to the correct Maxwell equations below. It appears in Eq. (30) because the standard Lagrangian is expressed in terms of $`E^2=E_x^2+E_y^2+E_z^2`$, whereas $`E_{\mathrm{}}^2`$ includes only the component of $`𝐄`$ in the direction of $`𝐝_{\mathrm{}}`$.

The electronic term in the Lagrangian, which includes the field-particle coupling, is

$`L_e`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}{\displaystyle \underset{i}{}}(c_i^{}\dot{c}_i\dot{c}_i^{}c_i){\displaystyle \underset{i,j}{}}c_i^{}H_{ij}c_j`$ (31a)
$`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}{\displaystyle \underset{i}{}}(b_i^{}\dot{b}_i\dot{b}_i^{}b_i){\displaystyle \underset{i,j}{}}b_i^{}\stackrel{~}{H}_{ij}b_j.`$ (31b)
Here $`H_{ij}=𝐱_i|H|𝐱_j`$ is the Hamiltonian in the absence of electromagnetic fields, while

$$\stackrel{~}{H}_{ij}=H_{ij}e^{i(\beta _i\beta _j)}+\mathrm{}\dot{\beta }_i\delta _{ij}.$$
(32)
The first expression (31a) for $`L_e`$ has exactly the same form as the Lagrangian in the case of no electromagnetic fields. This expresses the fundamental physical content of the principle of local gauge symmetry—that the influence of the field upon the particle can be expressed entirely in terms of the nonintegrable phase of the probability amplitude $`c_i=b_ie^{i\beta _i}`$. In the second expression (31b) for $`L_e`$, all of the nonintegrable phases are collected together in the effective Hamiltonian $`\stackrel{~}{H}_{ij}`$. This is the usual approach, in which the probability amplitude $`b_i`$ has a well-defined phase, and the field appears only in the Hamiltonian.

The Hamiltonian (32) appearing in the Lagrangian (31b) depends upon the phase difference $`\beta _i\beta _j`$. This phase difference is not well defined unless the sites $`i`$ and $`j`$ are linked. But according to postulate (i) above, all physical predictions of the theory must be unambiguous. Hence, the principle of local gauge symmetry demands that

$$H_{ij}=0(i\text{ not linked to }j).$$
(33)
Local gauge symmetry therefore imposes constraints not found in conventional tight-binding models.

Note that the Lagrangian $`L`$ is gauge invariant by construction. In other words, both $`L_f`$ and $`L_e`$ are invariant under the gauge (phase) transformation

$`b_i`$ $``$ $`b_ie^{i\chi _i},`$
$`\lambda _i`$ $``$ $`\lambda _i+\dot{\chi }_i,`$ (34)
$`\kappa _{ij}`$ $``$ $`\kappa _{ij}+\chi _i\chi _j,`$
where $`\chi _i`$ is an arbitrary integrable function.

Given the above Lagrangian, the Euler-Lagrange equation for $`\lambda _i`$ or $`\varphi _i`$ is just Gauss’s law

$$\underset{j}{}E_{ji}S_{ji}=4\pi q_i,$$
(35)
where $`q_i=eb_i^{}b_i`$ is the charge on site $`i`$. The corresponding equation for $`\kappa _{\mathrm{}}`$ is the Ampère-Maxwell equation

$$\underset{d_qS_{\mathrm{}}}{}B_qd_q\frac{1}{c}\frac{\text{d}}{\text{d}t}(E_{\mathrm{}}S_{\mathrm{}})=\frac{4\pi }{c}I_{\mathrm{}},$$
(36)
where $`I_{\mathrm{}}=I_{ji}=(2e/\mathrm{})\text{Im}(b_j^{}\stackrel{~}{H}_{ji}b_i)`$ is the current from site $`i`$ to site $`j`$. Summing (36) over all links that contain a given site $`i`$ yields the charge conservation law not (c)

$$\dot{q}_i+\underset{j}{}I_{ji}=0.$$
(37)
Thus, we see that $`\lambda _i`$ and $`\kappa _{\mathrm{}}`$ can be given a consistent interpretation as discrete electromagnetic potentials, since the above equations are in full agreement with macroscopic (i.e., long-wavelength) electromagnetism.

In some applications of Voronoi polyhedra, it is useful to link sites $`i`$ and $`j`$ whose polyhedra share only a line or point, hence $`S_{\mathrm{}}=S_{ij}=0`$. Okabe et al. (2000) For such links, the link volume $`\mathrm{\Omega }_{\mathrm{}}=\frac{1}{3}S_{\mathrm{}}d_{\mathrm{}}`$ is zero, so the electric field $`E_{\mathrm{}}`$ does not contribute to the field Lagrangian (30), Gauss’s law (35), or the Ampère-Maxwell equation (36). The magnetic-field contribution to (36) likewise vanishes, because $`S_{\mathrm{}}=0`$. The current through such a link must therefore be zero, which can only be true in general if the Hamiltonian matrix element $`H_{ij}`$ vanishes. Links with $`S_{ij}=0`$ are consequently devoid of any physical significance, and there is no loss of generality if one excludes them at the outset by linking only sites with $`S_{ij}>0`$.

Returning to the Lagrangian $`L`$, the Euler-Lagrange equation for $`b_i^{}`$ is just the Schrödinger equation

$$i\mathrm{}\dot{b}_i=\underset{j}{}\stackrel{~}{H}_{ij}b_j.$$
(38)
Since $`b_i`$ is an ordinary probability amplitude with a well-defined phase, $`\stackrel{~}{H}_{ij}`$ must be the Hamiltonian in the presence of electromagnetic fields. With the restriction (33), one can express (32) as

$`\stackrel{~}{H}_{ij}`$ $`=`$ $`H_{ij}\mathrm{exp}(i\kappa _{ij})+\mathrm{}\lambda _i\delta _{ij}`$ (39)
$`=`$ $`H_{ij}\mathrm{exp}(ieA_{ij}d_{ij}/\mathrm{}c)+e\varphi _i\delta _{ij}.`$
Note the strong similarity between this result and the Peierls substitution (4). The main difference is that (39) gives the Hamiltonian in the coordinate representation, not the tight-binding representation.

If $`|\kappa _{ij}|1`$ (i.e., if the field is weak or the lattice spacing is small), the Hamiltonian (39) reduces to

$$\stackrel{~}{H}_{ij}H_{ij}\frac{e}{mc}𝐀_{ij}𝐩_{ij}+\frac{e^2A_{ij}^2}{2mc^2}\mathrm{\Delta }_{ij}+e\varphi _i\delta _{ij}.$$
(40)
Here a vector potential has been defined as $`𝐀_{ij}=A_{ij}\widehat{𝐝}_{ij}`$, while the momentum operator is given by

$$𝐩_{ij}=\frac{m}{i\mathrm{}}𝐝_{ij}H_{ij}=\frac{m}{i\mathrm{}}(𝐱_i𝐱_j)H_{ij},$$
(41)
which is the same as the kinematic momentum defined above in Eq. (1). The Hamiltonian (39) therefore clearly gives the correct first-order $`𝐀𝐩`$ coupling. We shall see below that the dimensionless quantity

$$\mathrm{\Delta }_{ij}=\frac{1}{i\mathrm{}}𝐝_{ij}𝐩_{ij}=\frac{m}{\mathrm{}^2}d_{ij}^2H_{ij}$$
(42)
can be viewed as a geometric weight factor that gives the correct $`A^2`$ coupling also.

### IV.2 Geometric definition of momentum and kinetic energy

Up to this point, little has been said about the structure of the Hamiltonian $`H_{ij}`$. Within the bounds of the restriction (33), $`H_{ij}`$ may be treated as an arbitrary fitting parameter. However, in some circumstances it may be desirable to reduce the number of fitting parameters by using a theoretical formula for $`H_{ij}`$ that would reproduce the Schrödinger equation in the limit of zero lattice spacing.

Let us start by considering the momentum operator, which will be defined in this section as the canonical momentum $`𝐩=i\mathrm{}`$. A discrete expression for the gradient operator may be obtained from the integral definition of the gradient: Arfken (1985)

$$f(𝐱)=\underset{\mathrm{\Omega }0}{lim}\left[\frac{1}{\mathrm{\Omega }}_\mathrm{\Omega }f(𝐱)\text{d}𝐒\right],$$
(43)
where $`\text{d}𝐒`$ is a surface element pointing outward from $`\mathrm{\Omega }`$. Now the limiting volume in a discrete lattice is the volume $`\mathrm{\Omega }_i`$ of the Voronoi polyhedron for site $`𝐱_i`$. On the surface $`S_{ij}`$ shared by $`\mathrm{\Omega }_i`$ and $`\mathrm{\Omega }_j`$, the value of $`f(𝐱)`$ may be taken to be $`\frac{1}{2}[f(𝐱_i)+f(𝐱_j)]`$. Hence, the discrete gradient may be defined as

$$f(𝐱_i)=\frac{1}{\mathrm{\Omega }_i}\underset{j}{}\frac{1}{2}[f(𝐱_i)+f(𝐱_j)]𝐒_{ji},$$
(44)
where $`𝐒_{ji}=S_{ji}\widehat{𝐝}_{ji}`$. Now since

$$\underset{j}{}𝐒_{ji}=_{\mathrm{\Omega }_i}\text{d}𝐒=0,$$
(45)
the term involving $`f(𝐱_i)`$ drops out, leaving only

$$f(𝐱_i)=\frac{1}{2\mathrm{\Omega }_i}\underset{j}{}f(𝐱_j)𝐒_{ji}.$$
(46)
An alternative derivation of this result is given in Eq. (17) of Ref. Christ et al., 1982c.

The canonical momentum operator $`𝐩`$ may therefore be defined as Sakurai (1994)

$`𝐗_i|𝐩|\phi `$ $`=`$ $`i\mathrm{}𝐗_i|\phi `$ (47a)
$`=`$ $`{\displaystyle \frac{i\mathrm{}}{2\mathrm{\Omega }_i}}{\displaystyle \underset{j}{}}𝐗_j|\phi 𝐒_{ji}.`$ (47b)
Here the basis kets $`|𝐗_i=\mathrm{\Omega }_i^{1/2}|𝐱_i`$ are chosen to satisfy “$`\delta `$-function” normalization

$$𝐗_i|𝐗_j=\frac{\delta _{ij}}{\mathrm{\Omega }_i},$$
(48)
in contrast to the usual kets $`|𝐱_i`$, which are normalized to unity \[see Eq. (19)\]. The normalization (48) is used here because it agrees (in the limit $`\mathrm{\Omega }_i0`$) with the $`\delta `$-function normalization of continuous coordinate eigenkets, upon which the definition (47a) is based. Sakurai (1994)

Substituting $`|\phi =|𝐗_j`$ in Eq. (47b) then gives

$$𝐗_i|𝐩|𝐗_j=\frac{i\mathrm{}𝐒_{ij}}{2\mathrm{\Omega }_i\mathrm{\Omega }_j},$$
(49)
which clearly satisfies

$$\underset{j}{}𝐗_i|𝐩|𝐗_j\mathrm{\Omega }_j=0.$$
(50)
Replacing $`|𝐗_i=\mathrm{\Omega }_i^{1/2}|𝐱_i`$ in (49) then yields the desired result

$$𝐱_i|𝐩|𝐱_j=\frac{i\mathrm{}𝐒_{ij}}{2\sqrt{\mathrm{\Omega }_i\mathrm{\Omega }_j}}.$$
(51)
Note that this matrix is Hermitian, because $`𝐒_{ji}=𝐒_{ij}`$. If the kets $`|𝐱_i`$ are used in Eq. (47a) above, a non-Hermitian canonical momentum is obtained.

There is some question as to whether this definition of $`𝐩`$ should be referred to as “canonical,” because it does not satisfy the canonical commutation relations. In a continuous coordinate basis, the canonical momentum satisfies

$$𝐱^{}|[x^\alpha ,p^\beta ]|𝐱^{\prime \prime }=i\mathrm{}\delta _{\alpha \beta }\delta (𝐱^{}𝐱^{\prime \prime }),$$
(52)
where $`\alpha `$ and $`\beta `$ are Cartesian components of the given vectors. The corresponding equation in the discrete basis is

$$𝐗_i|[x^\alpha ,p^\beta ]|𝐗_j=\frac{i\mathrm{}d_{ij}^\alpha S_{ij}^\beta }{2\mathrm{\Omega }_i\mathrm{\Omega }_j},$$
(53)
which obviously does not agree. Note, however, that

$`{\displaystyle \underset{i,j}{}}𝐗_i|[x^\alpha ,p^\beta ]|𝐗_j\mathrm{\Omega }_i\mathrm{\Omega }_j`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}{\displaystyle \underset{i,j}{}}d_{ij}^\alpha S_{ij}^\beta `$ (54)
$`=`$ $`i\mathrm{}\delta _{\alpha \beta }\mathrm{\Omega },`$
where $`\mathrm{\Omega }`$ is the total volume, and the second equality is proved in Eq. (11) of Ref. Christ et al., 1982c. This agrees with the relation

$$𝐱^{}|[x^\alpha ,p^\beta ]|𝐱^{\prime \prime }\text{d}^3x^{}\text{d}^3x^{\prime \prime }=i\mathrm{}\delta _{\alpha \beta }\mathrm{\Omega }$$
(55)
in the continuous basis. Hence, Eq. (54) is as close as one can come to a canonical commutation relation in a general discrete basis. not (d, e)

A similar definition may be used for the kinetic energy operator $`T=\mathrm{}^2^2/2m`$. The integral definition of the divergence, Arfken (1985)

$$𝐅(𝐱)=\underset{\mathrm{\Omega }0}{lim}\left[\frac{1}{\mathrm{\Omega }}_\mathrm{\Omega }𝐅(𝐱)\text{d}𝐒\right],$$
(56)
gives the Laplacian

$$^2f(𝐱)=\underset{\mathrm{\Omega }0}{lim}\left[\frac{1}{\mathrm{\Omega }}_\mathrm{\Omega }f(𝐱)\text{d}𝐒\right],$$
(57)
the discrete form of which is

$$^2f(𝐱_i)=\frac{1}{\mathrm{\Omega }_i}\underset{j}{}\left(\frac{f(𝐱_j)f(𝐱_i)}{d_{ji}}\right)S_{ji}.$$
(58)
An alternative derivation of this result is given in Eq. (12) of Ref. Christ et al., 1982c. The procedure used above in Eqs. (47)–(51) then yields the kinetic energy operator

$$𝐱_i|T|𝐱_j=\frac{\mathrm{}^2}{2m}\frac{S_{ij}}{d_{ij}\sqrt{\mathrm{\Omega }_i\mathrm{\Omega }_j}}+\delta _{ij}\frac{\mathrm{}^2}{2m\mathrm{\Omega }_i}\underset{k}{}\frac{S_{ik}}{d_{ik}},$$
(59)
which satisfies \[cf. Eq. (50)\]

$$\underset{j}{}𝐗_i|T|𝐗_j\mathrm{\Omega }_j=0.$$
(60)
Note that for $`ij`$, $`𝐱_i|T|𝐱_j`$ decreases continuously to zero when $`S_{ij}0`$. This ensures that the Hamiltonian is a continuous function of the lattice coordinates, even as new links are formed and old ones are broken.

Such continuity is also desirable when the Hamiltonian is determined empirically, especially for applications (such as molecular dynamics Wang and Ho (1996); Xu et al. (1992)) in which the atomic positions vary with time. This can be achieved by defining the nonlocal elements of the empirical Hamiltonian as

$$𝐱_i|H|𝐱_j=\frac{\mathrm{}^2}{2m}\frac{S_{ij}}{d_{ij}\sqrt{\mathrm{\Omega }_i\mathrm{\Omega }_j}}f_{ij}(ij),$$
(61)
where the fitting parameter $`f_{ij}`$ is a continuous, nonsingular function of the lattice coordinates.

Note that the operators $`𝐩`$ and $`T`$ do not satisfy $`T=𝐩^2/2m`$, because $`𝐩^2`$, unlike $`𝐩`$ and $`T`$, couples sites that are not linked. However, $`𝐩`$ and $`T`$ are related by

$$𝐩_{ij}=\frac{m}{i\mathrm{}}𝐝_{ij}T_{ij}$$
(62a)
or
$$𝐩=\frac{m}{i\mathrm{}}[𝐱,T].$$
(62b)
Thus, any Hamiltonian of the form $`H=T+V(𝐱)`$, where $`V`$ is a local potential, satisfies Eq. (1). Hence, for such a Hamiltonian, the canonical momentum $`𝐩=i\mathrm{}`$ used in this section agrees with the kinematic momentum defined earlier.

Now let us examine the dimensionless factor $`\mathrm{\Delta }_{ij}`$ defined above in Eq. (42). If $`H=T+V`$, this becomes

$$\mathrm{\Delta }_{ij}=\frac{S_{ij}d_{ij}}{2\sqrt{\mathrm{\Omega }_i\mathrm{\Omega }_j}}=\frac{3\mathrm{\Omega }_{ij}}{2\sqrt{\mathrm{\Omega }_i\mathrm{\Omega }_j}},$$
(63)
where $`\mathrm{\Omega }_{ij}=\frac{1}{3}S_{ij}d_{ij}`$ is the volume of the link between sites $`i`$ and $`j`$ (see Appendix B). The factor $`\mathrm{\Delta }_{ij}`$ appears in the $`A^2`$ term in the Hamiltonian (40), which will be referred to as $`H_2`$. In a continuous coordinate basis, $`H_2`$ is given by

$$𝐱^{}|H_2|𝐱^{\prime \prime }=\frac{e^2}{2mc^2}A^2(𝐱^{})\delta (𝐱^{}𝐱^{\prime \prime }),$$
(64)
which means that it satisfies

$$𝐱^{}|H_2|𝐱^{\prime \prime }\text{d}^3x^{}\text{d}^3x^{\prime \prime }=\frac{e^2}{2mc^2}A^2(𝐱^{})\text{d}^3x^{}.$$
(65)
The corresponding equations for the discrete basis are

$$𝐗_i|H_2|𝐗_j=\frac{e^2}{2mc^2}\frac{3A_{ij}^2\mathrm{\Omega }_{ij}}{2\mathrm{\Omega }_i\mathrm{\Omega }_j}$$
(66)
and

$`{\displaystyle \underset{i,j}{}}𝐗_i|H_2|𝐗_j\mathrm{\Omega }_i\mathrm{\Omega }_j`$ $`=`$ $`{\displaystyle \frac{e^2}{2mc^2}}{\displaystyle \frac{3}{2}}{\displaystyle \underset{i,j}{}}A_{ij}^2\mathrm{\Omega }_{ij}`$ (67)
$`=`$ $`{\displaystyle \frac{e^2}{2mc^2}}{\displaystyle \underset{\mathrm{}}{}}3A_{\mathrm{}}^2\mathrm{\Omega }_{\mathrm{}}.`$
The second equality in (67) was obtained by noting that a sum over $`i`$ and $`j`$ covers each link $`\mathrm{}`$ twice. The only apparent difference between Eqs. (65) and (67) is a factor of 3. This appears for the same reason that it does in the Lagrangian (30)—i.e., $`A^2`$ in Eq. (65) refers to $`A_x^2+A_y^2+A_z^2`$, whereas $`A_{\mathrm{}}^2`$ in Eq. (67) refers only to the component of $`𝐀`$ in the direction of $`𝐝_{\mathrm{}}`$.

Therefore, Eqs. (65) and (67) are the same in the limit of zero lattice spacing, and the factor $`\mathrm{\Delta }_{ij}`$ is simply a geometric weight factor that provides the correct $`A^2`$ coupling in the Hamiltonian (40).

### IV.3 Spin

The theory presented thus far has been for a particle with spin zero. Particles with spin $`\frac{1}{2}`$ may be described using a discrete version of the Dirac Hamiltonian for a free particle:

$$H=c𝜶𝐩+\beta mc^2,$$
(68)
where $`𝜶`$ and $`\beta `$ are Dirac’s $`4\times 4`$ matrices. The momentum operator $`𝐩`$ can either be calculated from geometry or fitted to experiment. In the presence of electromagnetic fields, the Hamiltonian becomes

$$\stackrel{~}{H}=c𝜶𝝅+\beta mc^2+e\varphi ,$$
(69)
where \[cf. Eq. (39)\]

$$𝝅_{ij}=𝐩_{ij}\mathrm{exp}(i\kappa _{ij}).$$
(70)
A nonrelativistic Hamiltonian may be obtained by applying a Foldy-Wouthuysen transformation Foldy and Wouthuysen (1950); Messiah (1962) to Eq. (69), which yields

$`H_{\text{nr}}`$ $`=`$ $`{\displaystyle \frac{1}{2m}}(𝝈𝝅)^2{\displaystyle \frac{1}{8m^3c^2}}(𝝈𝝅)^4+e\varphi `$
$`{\displaystyle \frac{1}{8m^2c^2}}[𝝈𝝅,([𝝈𝝅,e\varphi ]+i\mathrm{}𝝈\dot{𝝅})],`$
where $`𝝈`$ is the Pauli spin matrix, and all terms of order $`(v/c)^4`$ have been included. This Hamiltonian couples sites that are not linked, but there is no ambiguity because the Dirac equation is taken as fundamental.

If we assume for simplicity that the lattice coordinates do not depend on time, then

$$([𝝈𝝅,e\varphi ]+i\mathrm{}𝝈\dot{𝝅})_{ij}=𝝈𝝅_{ij}V_{ij},$$
(72)
where

$$V_{ij}=e(\varphi _i\varphi _j)\mathrm{}\dot{\kappa }_{ij}=eE_{ij}d_{ij}$$
(73)
is the difference in potential energy of sites $`i`$ and $`j`$ due to the electric field. The last term in Eq. (IV.3) therefore consists of the Darwin term

$$H_{ij}^\text{D}=\frac{1}{8m^2c^2}\underset{k}{}(V_{ki}+V_{kj})𝝅_{ik}𝝅_{kj}$$
(74)
plus the spin-orbit coupling

$$H_{ij}^{\text{so}}=\frac{i}{8m^2c^2}\underset{k}{}(V_{ki}+V_{kj})𝝈(𝝅_{ik}\times 𝝅_{kj}),$$
(75)
where the identity

$$(𝝈𝐚)(𝝈𝐛)=𝐚𝐛+i𝝈(𝐚\times 𝐛)$$
(76)
has been used. Now the main contribution to spin-orbit coupling comes from the atomic cores, where the potential energy and wave function vary rapidly. However, in any basis of reasonable size, the lattice imposes a wavelength cutoff that eliminates such rapid variations. The potential $`\varphi _i`$ must therefore be viewed as a pseudopotential, not a true atomic potential. Hence, for practical purposes, $`\varphi _i`$ in the spin-orbit Hamiltonian (75) should be treated as a fitting parameter that is independent of the value used for the first $`e\varphi `$ term in (IV.3).

The first two terms in the Hamiltonian (IV.3) are kinetic-energy terms, which may be rewritten using

$$(𝝈𝝅)_{ij}^2=\underset{k}{}[𝝅_{ik}𝝅_{kj}+i𝝈(𝝅_{ik}\times 𝝅_{kj})],$$
(77)
in which the second term describes the intrinsic magnetic dipole moment of the particle. For a general lattice, this term is not zero even when there is no electromagnetic field, because different components of the momentum operator do not commute (i.e., $`𝐩\times 𝐩0`$). This follows from the fact that there is generally no more than one site $`k`$ linked to both $`i`$ and $`j`$, and for that $`i`$ and $`j`$, $`𝐩_{ik}\times 𝐩_{kj}`$ is generally not zero.

However, if the lattice is a Bravais lattice, then $`𝐩\times 𝐩`$ is always zero. This follows from the fact that every site in a Bravais lattice is identical, so for a given nonzero $`𝐩_{ik}\times 𝐩_{kj}`$, there is always another site $`l`$ with $`𝐩_{il}=𝐩_{kj}`$ and $`𝐩_{lj}=𝐩_{ik}`$, hence $`𝐩_{ik}\times 𝐩_{kj}+𝐩_{il}\times 𝐩_{lj}=0`$. Clearly one also has $`𝐝_{ik}\times 𝐝_{kj}+𝐝_{il}\times 𝐝_{lj}=0`$, so the sites $`i`$, $`k`$, $`j`$, and $`l`$ lie in a single plane. If $`i`$ is not linked to $`j`$ and $`k`$ is not linked to $`l`$, then $`i`$, $`k`$, $`j`$, and $`l`$ form a single plaquette $`q`$, which has the shape of a rectangle. Otherwise, they form two triangular plaquettes.

If the momentum operator is given by Eq. (51), then the intrinsic magnetic dipole term in the Hamiltonian (IV.3) for such a Bravais lattice is

$$H_{ij}^{\text{mag}}=\frac{e\mathrm{}}{8mc}\left(\frac{S_{ik}S_{kj}S_q^2}{d_{ik}d_{kj}\mathrm{\Omega }_i^2}\right)𝝈𝐁_q,$$
(78)
where the weak-field approximation $`|\kappa _{ij}|1`$ has been used, and the direction of $`𝐁_q`$ is that of $`𝐒_q`$. If the sites $`i`$, $`k`$, $`j`$, and $`l`$ form a single rectangular plaquette, then $`S_q`$ is the area of that plaquette; otherwise, it is the combined area of the two triangles (in which case $`𝐁_q`$ is the average magnetic field of the two plaquettes).

As an example, consider a simple cubic lattice with lattice constant $`a`$, for which $`S_{ik}=S_{kj}=S_q=a^2`$, $`d_{ik}=d_{kj}=a`$, and $`\mathrm{\Omega }_i=a^3`$. In this case, the factor in parentheses in Eq. (78) is unity, and $`H_{ij}^{\text{mag}}`$ couples sites on opposite corners of each plaquette (with $`d_{ij}=\sqrt{2}a`$). By comparison, the dipole term in the continuum Hamiltonian is given by

$$\frac{1}{2m}[𝝈(𝐩e𝐀/c)]^2=\frac{1}{2m}(𝐩e𝐀/c)^2\frac{e\mathrm{}}{2mc}𝝈𝐁.$$
(79)
The numerical factor in front of this dipole coupling is four times larger than that in Eq. (78). This occurs because (78) couples each site $`i`$ to four other sites $`j`$.

## V Application to tetrahedral semiconductors

This section considers several different methods of implementing the theory developed in Sec. IV. Spin is neglected in all of the applications that follow.

### V.1 Discrete pseudopotential method

The simplest geometry occurs when the lattice sites $`𝐱_i`$ are chosen to lie on a Bravais lattice. One possible approach in this case is to use the geometric expression (59) for the kinetic energy $`T`$, and assume that the potential energy $`V`$ is local. This approach will be referred to as the discrete pseudopotential method.

If $`𝐱`$ lies on a Bravais lattice (the subscript $`i`$ is omitted here), one may define a reciprocal lattice as the set of all vectors $`𝐠`$ such that $`𝐠𝐱=2\pi \times \text{integer}`$. The volume of a primitive cell in the direct lattice is denoted $`\omega _0`$, while that of a primitive cell in the reciprocal lattice is $`\omega _0^{}=(2\pi )^3/\omega _0`$.

In a crystal, the Hamiltonian will be periodic with respect to some larger Bravais lattice whose sites are denoted by $`𝐑`$, where $`𝐑\{𝐱\}`$. One may then define a corresponding reciprocal lattice as the set of all vectors $`𝐆`$ such that $`𝐆𝐑=2\pi \times \text{integer}`$. The volume of a primitive cell in $`𝐑`$ space is $`\mathrm{\Omega }_0=n\omega _0`$, where $`n`$ is an integer, and the volume of a primitive cell in $`𝐆`$ space is $`\mathrm{\Omega }_0^{}=(2\pi )^3/\mathrm{\Omega }_0`$.

Periodic boundary conditions may then be implemented over an even larger Bravais lattice whose sites are denoted by $`𝐋`$, where $`𝐋\{𝐑\}`$. The corresponding reciprocal lattice vectors are denoted $`𝐤`$, where $`𝐤𝐋=2\pi \times \text{integer}`$. The volume over which periodic boundary conditions are applied is $`\mathrm{\Omega }=N\mathrm{\Omega }_0`$, where $`N`$ is an integer, and the volume of a primitive cell in $`𝐤`$ space is $`\mathrm{\Omega }^{}=(2\pi )^3/\mathrm{\Omega }`$. Note that according to the above definitions, $`𝐆\{𝐤\}`$ and $`𝐠\{𝐆\}`$.

Coordinate eigenkets in this system are denoted $`|𝐱`$, where $`|𝐱+𝐋|𝐱`$ due to the periodic boundary conditions. The orthogonality and closure relations in this basis are therefore

$`𝐱|𝐱^{}={\displaystyle \underset{𝐋}{}}\delta _{𝐱𝐱^{},𝐋},`$ (80)
$`{\displaystyle \underset{𝐱\mathrm{\Omega }}{}}|𝐱𝐱|=1.`$ (81)
In a system with periodic boundary conditions, the coordinate operator is not well defined; only periodic functions of the coordinate are permitted. Therefore, the definition of the kinematic momentum must be slightly modified. Instead of (1), one has

$$𝐱|𝐩|𝐱^{}=\frac{m}{i\mathrm{}}(𝐱𝐱^{})𝐱|H|𝐱^{}\text{if }𝐱\mathrm{\Omega }\text{ and }𝐱^{}\mathrm{\Omega },$$
(82)
with $`𝐱|𝐩|𝐱^{}=𝐱+𝐋|𝐩|𝐱^{}`$ otherwise.

Another useful representation is the crystal momentum representation

$$|𝐤=\frac{1}{\sqrt{𝒩}}\underset{𝐱\mathrm{\Omega }}{}e^{i𝐤𝐱}|𝐱,$$
(83)
where $`𝒩=nN=\mathrm{\Omega }/\omega _0`$. The corresponding orthogonality and closure relations are

$`𝐤|𝐤^{}={\displaystyle \underset{𝐠}{}}\delta _{𝐤𝐤^{},𝐠},`$ (84)
$`{\displaystyle \underset{𝐤\omega _0^{}}{}}|𝐤𝐤|=1.`$ (85)
In the crystal momentum representation, a periodic Hamiltonian $`𝐱|H|𝐱^{}=𝐱+𝐑|H|𝐱^{}+𝐑`$ couples only those states that differ by a reciprocal lattice vector $`𝐆`$:

$$𝐤^{}|H|𝐤=\underset{𝐆}{}\delta _{𝐤^{},𝐤+𝐆}𝐤+𝐆|H|𝐤,$$
(86)
where

$$𝐤+𝐆|H|𝐤=\frac{1}{n}\underset{𝐱\mathrm{\Omega }}{}\underset{𝐱^{}\mathrm{\Omega }_0}{}e^{i(𝐤+𝐆)𝐱^{}}𝐱^{}|H|𝐱e^{i𝐤𝐱}.$$
(87)
The kinematic momentum (82) also satisfies Eq. (86). Its matrix elements that are related to those of $`H`$ by

$$𝐤+𝐆|𝐩|𝐤=\frac{m}{\mathrm{}}_𝐤𝐤+𝐆|H𝐤.$$
(88)
On a Bravais lattice, the kinetic energy (59) and canonical momentum (51) are translationally invariant:

$$𝐱|T|𝐱^{}=T(𝐱𝐱^{}),$$
(89)
where $`T(𝐱+𝐋)=T(𝐱)`$. The matrix elements of $`T`$ are therefore given by

$$𝐤|T|𝐤^{}=T(𝐤)\underset{𝐠}{}\delta _{𝐤𝐤^{},𝐠},$$
(90)
where

$$T(𝐤)=\underset{𝐱\mathrm{\Omega }}{}T(𝐱)e^{i𝐤𝐱}.$$
(91)
For tetrahedral semiconductors with the diamond or zinc-blende structure, it is convenient to use a cubic lattice for the grid $`\{𝐱\}`$. Expressions for the link distances $`d_{ij}`$ and surface areas $`S_{ij}`$ are given in Appendix C for the simple cubic, body-centered cubic, and face-centered cubic lattices. The resulting kinetic-energy operators given by Eqs. (59) and (91) are

$`T_{\text{sc}}(𝐤)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{}^2}{ma^2}}[\mathrm{sin}^2(\frac{1}{2}k_xa)+\mathrm{sin}^2(\frac{1}{2}k_ya)+\mathrm{sin}^2(\frac{1}{2}k_za)],`$
$`T_{\text{bcc}}(𝐤)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2ma^2}}\{6[1\mathrm{cos}(\frac{1}{2}k_xa)\mathrm{cos}(\frac{1}{2}k_ya)\mathrm{cos}(\frac{1}{2}k_za)]`$
$`+[\mathrm{sin}^2(\frac{1}{2}k_xa)+\mathrm{sin}^2(\frac{1}{2}k_ya)+\mathrm{sin}^2(\frac{1}{2}k_za)]\},`$
$`T_{\text{fcc}}(𝐤)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{}^2}{ma^2}}[3\mathrm{cos}(\frac{1}{2}k_ya)\mathrm{cos}(\frac{1}{2}k_za)`$ (92)
$`\mathrm{cos}(\frac{1}{2}k_za)\mathrm{cos}(\frac{1}{2}k_xa)`$
$`\mathrm{cos}(\frac{1}{2}k_xa)\mathrm{cos}(\frac{1}{2}k_ya)],`$
all of which reduce to $`T(𝐤)\mathrm{}^2k^2/2m`$ when $`ka1`$. Here $`a`$ is the lattice constant of the grid $`\{𝐱\}`$, which is some integer fraction of the lattice constant $`a_0`$ of the crystal lattice $`\{𝐑\}`$.

A canonical momentum operator $`𝐩(𝐤)`$ corresponding to Eq. (51) may be defined in a similar manner. This operator is given by

$$𝐩(𝐤)=\frac{m}{\mathrm{}}_𝐤T(𝐤),$$
(93)
which follows from Eq. (62). Note that this result is just a special case of the kinematic momentum (88).

The matrix elements of a local periodic potential $`V(𝐱)=V(𝐱+𝐑)`$ are given by Eq. (86), where

$$𝐤+𝐆|V|𝐤=\frac{1}{n}\underset{𝐱\mathrm{\Omega }_0}{}V(𝐱)e^{i𝐆𝐱}$$
(94)
is independent of $`𝐤`$. It will be assumed here that $`V`$ is a superposition of local atomic pseudopotentials:

$$V(𝐱)=\underset{\mu =1}{\overset{N_a}{}}\underset{𝐑\mathrm{\Omega }}{}v_\mu (𝐱𝐑𝝉_\mu ),$$
(95)
where $`v_\mu (𝐱)=v_\mu (𝐱+𝐋)`$ is the pseudopotential for atom $`\mu `$, whose position in the unit cell $`\mathrm{\Omega }_0`$ is given by $`𝝉_\mu `$. In this case

$$𝐤+𝐆|V|𝐤=\frac{1}{N_a}\underset{\mu =1}{\overset{N_a}{}}v_\mu (𝐆)e^{i𝐆𝝉_\mu },$$
(96)
where $`v_\mu (𝐆)`$ is the atomic form factor

$$v_\mu (𝐆)=\frac{N_a}{n}\underset{𝐱\mathrm{\Omega }}{}v_\mu (𝐱)e^{i𝐆𝐱},$$
(97)
and $`N_a`$ is the number of atoms in the unit cell $`\mathrm{\Omega }_0`$.

The main practical difficulty in implementing the discrete pseudopotential method is that $`T(𝐤)`$ is not a good approximation to the continuum kinetic energy

$$T_{\text{cont}}(𝐤)=\frac{\mathrm{}^2k^2}{2m}$$
(98)
unless $`ka1`$. This means that the lattice constant $`a`$ of the grid $`\{𝐱\}`$ must be significantly smaller than the lattice constant $`a_0`$ of the crystal lattice $`\{𝐑\}`$. To obtain one grid point at each atom in the diamond structure, $`a`$ must satisfy $`a=a_0/2l`$ (for a bcc grid) or $`a=a_0/4l`$ (for sc and fcc), where $`l`$ is a positive integer. Numerical accuracy generally requires $`l>1`$, as shown below.

The shape of the Brillouin zone is also important. If the shape of the Wigner-Seitz cell for $`\omega _0^{}`$ is not congruent with the shape of the Wigner-Seitz cell for $`\mathrm{\Omega }_0^{}`$, then $`T(𝐤)`$ will deviate from $`T_{\text{cont}}(𝐤)`$ more rapidly in some directions than others. This can lead to significant qualitative errors in the kinetic energy. For example, in diamond, the lowest continuum eigenvalues at $`\mathrm{\Gamma }`$ are $`T_{\text{cont}}(𝐆)=\mathrm{}^2G^2/2m`$, where $`𝐆=(000)`$, $`111`$, $`200`$, and $`220`$ (in units of $`2\pi /a_0`$). But for a sc grid with $`a=\frac{1}{4}a_0`$, the value of $`T_{\text{sc}}(200)=2\mathrm{}^2/ma^2`$ is actually lower than that of $`T_{\text{sc}}(111)=3\mathrm{}^2/ma^2`$. The correct ratio $`T_{\text{cont}}(200)=\frac{4}{3}T_{\text{cont}}(111)`$ is only approached in the limit of very small $`a/a_0`$, making the sc grid a poor choice for diamond.

The natural choice for diamond is the fcc grid, since its Brillouin zone has the same shape as that of diamond. Indeed, one has $`T_{\text{fcc}}(200)=\frac{4}{3}T_{\text{fcc}}(111)`$ even for the maximum grid size $`a=\frac{1}{4}a_0`$. The only problem here is that $`T_{\text{fcc}}(220)=\frac{3}{2}T_{\text{fcc}}(200)`$, which is not sufficiently close to the correct ratio $`T_{\text{cont}}(220)=2T_{\text{cont}}(200)`$ to make $`a=\frac{1}{4}a_0`$ a satisfactory choice for numerical calculations. The next possibility is $`a=\frac{1}{8}a_0`$, which yields $`T_{\text{fcc}}(220)=(\frac{3}{2}+\frac{\sqrt{2}}{4})T_{\text{fcc}}(200)1.85T_{\text{fcc}}(200)`$.

The energy band structure for GaAs calculated using an fcc grid with $`a=\frac{1}{8}a_0`$ is given in Fig. 1.

The fcc kinetic energy obtained from Eq. (92) was multiplied by a constant factor $`\pi ^2(2+\sqrt{2})/32=1.053`$ so that $`T_{\text{fcc}}(𝐤)`$ matches $`T_{\text{cont}}(𝐤)`$ at $`𝐆=111`$ and $`𝐆=200`$. The pseudopotential form factors (97) for this calculation were taken from Ref. Cohen and Bergstresser, 1966. No attempt was made to fit the energy bands by modifying the form factors; the purpose of this figure is merely to demonstrate the close similarity between the discrete fcc band structure and the continuum band structure. Slight adjustments in the model parameters would likely give an even better agreement.

The main problem with this result is that $`a=\frac{1}{8}a_0`$ corresponds to a basis size of 512 grid points per primitive unit cell $`\mathrm{\Omega }_0`$. This is unattractive in comparison to the basis dimensions of approximately 100 plane waves that are typically used in empirical pseudopotential calculations for tetrahedral semiconductors. However, changing the fcc grid to $`a=\frac{1}{4}a_0`$ (i.e., 64 grid points per primitive cell) makes it impossible to achieve a satisfactory fit to the band structure using local pseudopotentials. A good fit is possible only if the nonlocal Hamiltonian matrix elements are treated as fitting parameters. But in that case, one can reduce the basis dimensions even further by using the tight-binding approach.

### V.2 Tight-binding method

In the tight-binding approach, the grid points are no longer restricted to lie on a Bravais lattice, and all of the Hamiltonian matrix elements are treated as fitting parameters. In addition, it is assumed here that the model is constructed using the symmetrization procedure described in Sec. II, so that a distinct set of orthogonal orbitals is associated with each atom. The objective then is to find the smallest coordinate basis that provides a physically reasonable model of a given system.

The basis kets in the tight-binding approach will be written as $`|\alpha ,𝐑+𝝉_\mu `$, where $`𝐑`$ is a lattice vector for the Bravais lattice over which the Hamiltonian is periodic, and $`𝝉_\mu `$ is the position of atom $`\mu `$ within the primitive unit cell $`\mathrm{\Omega }_0`$. These quantities are defined exactly as they were in Sec. V.1; the vectors $`𝐋`$, $`𝐆`$, and $`𝐤`$ are also defined in the same way. The label $`\alpha `$ may refer to a coordinate $`𝐱_\alpha `$ within the atom, in which case

$$|\alpha ,𝐑+𝝉_\mu |𝐱_\alpha +𝐑+𝝉_\mu $$
(99)
is just a coordinate eigenket. However, $`\alpha `$ may also be used as a symmetry label for an atomic orbital that is a symmetrized linear combination of the kets (99):

$$|\alpha ,𝐑+𝝉_\mu =\underset{\beta }{}C_\beta (\alpha )|𝐱_\beta +𝐑+𝝉_\mu .$$
(100)
In either case, the basis is orthogonal:

$$\alpha ,𝐑+𝝉_\mu |\alpha ^{},𝐑^{}+𝝉_\mu ^{}=\delta _{\alpha \alpha ^{}}\delta _{\mu \mu ^{}}\underset{𝐋}{}\delta _{𝐑𝐑^{},𝐋},$$
(101)
and complete:

$$\underset{\alpha }{}\underset{\mu }{}\underset{𝐑\mathrm{\Omega }}{}|\alpha ,𝐑+𝝉_\mu \alpha ,𝐑+𝝉_\mu |=1.$$
(102)
In periodic systems, it is convenient to define the Bloch sums Bloch (1928)

$$|\alpha ,\mu ,𝐤=\frac{1}{\sqrt{N}}\underset{𝐑\mathrm{\Omega }}{}e^{i𝐤(𝐑+𝝉_\mu )}|\alpha ,𝐑+𝝉_\mu ,$$
(103)
which are also orthogonal and complete:

$$\alpha ,\mu ,𝐤|\alpha ^{},\mu ^{},𝐤^{}=\delta _{\alpha \alpha ^{}}\delta _{\mu \mu ^{}}\underset{𝐆}{}\delta _{𝐤𝐤^{},𝐆}e^{i𝐆𝝉_\mu },$$
(104)
$$\underset{\alpha }{}\underset{\mu }{}\underset{𝐤\mathrm{\Omega }_0^{}}{}|\alpha ,\mu ,𝐤\alpha ,\mu ,𝐤|=1.$$
(105)
If the Hamiltonian $`H`$ is invariant with respect to lattice translations $`𝐑`$, then its matrix elements in the Bloch basis are

$`\alpha ,\mu ,𝐤|H|\alpha ^{},\mu ^{},𝐤^{}`$ $`=`$ $`\alpha ,\mu ,𝐤|H|\alpha ^{},\mu ^{},𝐤`$
$`\times {\displaystyle \underset{𝐆}{}}\delta _{𝐤𝐤^{},𝐆}e^{i𝐆𝝉_\mu ^{}},`$
where

$`\alpha ,\mu ,𝐤|H|\alpha ^{},\mu ^{},𝐤`$ $`=`$ $`{\displaystyle \underset{𝐑^{}\mathrm{\Omega }}{}}\alpha ,𝝉_\mu |H|\alpha ^{},𝐑^{}+𝝉_\mu ^{}`$ (107)
$`\times e^{i𝐤(𝐑^{}+𝝉_\mu ^{}𝝉_\mu )}.`$
The kinematic momentum operator (41) is related to this Hamiltonian by

$`\alpha ,\mu ,𝐤|𝐩|\alpha ^{},\mu ^{},𝐤`$ $`=`$ $`{\displaystyle \frac{m}{i\mathrm{}}}(𝐱_\alpha 𝐱_\alpha ^{}+i_𝐤)`$ (108)
$`\times \alpha ,\mu ,𝐤|H|\alpha ^{},\mu ^{},𝐤,`$
where the original basis was assumed to be given by (99). Note that (108) has the form of an intra-atomic matrix element (proportional to $`𝐱_\alpha 𝐱_\alpha ^{}`$) plus an inter-atomic matrix element (proportional to $`i_𝐤`$). The zero-parameter model of Eq. (3) is obtained in the limit $`𝐱_\alpha 0`$.

Let us now consider specific examples of $`H`$ for tetrahedral semiconductors. The simplest tight-binding model for the diamond or zinc-blende structure consists of a single $`s`$ orbital per atom, which is obtained by putting one coordinate eigenket at each atomic position \[see Eq. (17)\]. In this model, as shown in Appendix C, each atom is linked to 4 nearest neighbors and 12 second-nearest neighbors. More distant linkages do not exist because the Voronoi polyhedra for these atoms do not touch one another.

Such a simple model is, of course, unable to describe even the qualitative features of tetrahedral semiconductors. The simplest conventional tight-binding model that works in this case is the $`sp^3`$ model. Slater and Koster (1954); Harrison (1989, 1999); Tserbak et al. (1993) The basic features of this model are described in Table 1, which lists the basis size (number of orbitals per atom) and the number of independent Hamiltonian matrix elements for coupling between atoms in the diamond structure out to third nearest neighbors. Tserbak et al. (1993)

Table 1 also lists the properties of other tight-binding models used in the literature, such as $`sp^3s^{}`$, Vogl et al. (1983) $`sp^3d^2`$, Chang and Aspnes (1990) and $`sp^3d^5s^{}`$. Jancu et al. (1998) The number of parameters listed in this table is the number permitted by the symmetry of the model, which is not necessarily the same as that used in any specific implementation in the literature.

For comparison, the bottom half of Table 1 lists the properties of several tight-binding models constructed from symmetrized coordinate eigenkets. The number of free parameters for the models generated by $`|100`$, $`|111`$, and $`|110`$ may be deduced easily from the link geometry results presented in Appendix C. The corresponding numbers for the compound models (with more than one generator) were determined using the algorithm in Appendix B.

The most striking feature of Table 1 is the relative paucity of free parameters in the symmetrized coordinate approach, which occurs because of the restriction (33) imposed by local gauge symmetry. Several of the symmetrized coordinate models are direct analogs of conventional models (i.e., they have identical symmetry); for example, the $`|111`$ ($`T_d`$) model corresponds to $`sp^3`$, and the $`|100`$ model corresponds to $`sp^3d^2`$ \[see Eqs. (14) and (17)\]. However, the number of free parameters in conventional tight-binding theory grows steadily with distance, whereas in the present theory there is no coupling beyond second nearest neighbors.

This dearth of adjustable parameters means that the smallest basis sets do not provide a reliable model for the energy band structure. For example, in the $`|111`$ ($`T_d`$) model, the splitting of the bonding and antibonding $`s`$ states at the $`\mathrm{\Gamma }`$ point is the same as that of the $`p`$ states at $`\mathrm{\Gamma }`$ (and that of the $`p`$ states at $`X`$). In the $`|100`$ model, there is no coupling at all between $`p`$ orbitals on different atoms at the $`\mathrm{\Gamma }`$ point, so the splitting of bonding and antibonding states is zero. Such difficulties arise primarily because there is only one nearest-neighbor coupling parameter in these models.

To increase the number of adjustable parameters without an undue increase in basis size, one must deliberately search for models with the most complicated topology available. A good starting point is the $`|111`$ ($`O_h`$) basis, which already has 3 nearest-neighbor parameters. Combining this with a $`|100`$ basis raises that number to 4 or 5, depending on the relative values of the coordinates in the two generators. This 14-orbital model is a dramatic improvement over any of the smaller models; however, an extra $`|000`$ site adds substantial extra flexibility without much change in the basis size. Thus, the 15-orbital model generated by $`|000`$, $`|111`$, and $`|100`$ is probably the smallest basis capable of describing tetrahedral semiconductors accurately. In the language of conventional tight-binding theory, this would be referred to as an $`s^3p^6d^5f`$ model.

As shown in Table 1, this model has 17 free parameters (one of which is just the reference energy). Specific definitions for these parameters are given in Table 2, which also presents parameter values for Ge and Si obtained by fitting the band structure of the 15-orbital model to that given by the nonlocal empirical pseudopotentials of Chelikowsky and Cohen. Chelikowsky and Cohen (1973, 1974, 1976)

Local Hamiltonian matrix elements are labeled $`V`$, whereas nonlocal on-site, nearest-neighbor, and second-nearest-neighbor terms are denoted $`\alpha `$, $`\beta `$, and $`\gamma `$, respectively. The subscripts $`a`$, $`b`$, $`e`$, and $`f`$ refer to independent sites generated by $`|0,0,0`$, $`|r,r,r`$, $`|r,r,r`$, and $`|r^{},0,0`$, respectively. The labels $`a`$, $`e`$, and $`f`$ are the Wyckoff labels for sites of different symmetry in the diamond structure. Hahn (1992) The label $`b`$ is not correct Wyckoff notation; it actually refers to an independent $`e`$ site, but the notation $`b`$ was used here because these sites lie on the bonds between atoms.

The energy band structure calculated from the parameters in Table 2 is plotted in Fig. 2.

One can see that the 15-orbital model provides a good fit to the nonlocal pseudopotential bands from the bottom of the valence band to about 5 eV above the top of the valence band. Qualitative errors begin to occur near 9 eV at both the $`L`$ and $`X`$ points. For example, the $`X_1`$ level near 9 eV in both figures should occur above 12 eV. This discrepancy can be eliminated, but the author has not found any way of doing so without adversely affecting the quality of the overall fit.

It should be emphasized that the parameters in Table 2 are presented here merely as “proof of concept;” they are in no way intended as the final word on the subject, and the author would be surprised if a better set were not found in the future. The quantities included in the fitting procedure were the valence- and conduction-band energy levels at $`\mathrm{\Gamma }`$, $`X`$, $`L`$, and $`K`$. Effective masses and deformation potentials were not considered, and no attempt was made at ensuring transferability.

The main difficulty encountered during the fitting was the lack of any reliable method for establishing a sound starting point. Unlike the case for smaller tight-binding models, the present Hamiltonian has almost no simple analytical solutions (except for the $`\mathrm{\Gamma }_{12^{}}`$, $`\mathrm{\Gamma }_{12}`$, and $`X_2`$ states, which are relatively unimportant) that can be used to determine starting values. The formula (59) for the kinetic energy provides a set of “free-particle” parameters that is better than nothing, but in a 15-orbital basis, Eq. (59) is a poor approximation to the continuum kinetic energy (98). After several months and dozens of different schemes (which still sampled only an infinitesimal fraction of parameter space), the author was unable to find any method whose success could honestly be attributed to anything other than trial and error. Hence, the development of a robust fitting procedure remains an unsolved problem.

### V.3 Dielectric function

As a test of the field-particle coupling in the 15-orbital model, the imaginary part of the transverse dielectric tensor was calculated from the formula Graf and Vogl (1995); Bassani and Pastori Parravicini (1975)

$`ϵ_2^{\alpha \beta }(\omega )`$ $`=`$ $`{\displaystyle \frac{4\pi ^2e^2}{m^2\omega ^2}}{\displaystyle \underset{c,v}{}}{\displaystyle \frac{2}{(2\pi )^3}}{\displaystyle _{\mathrm{\Omega }_0^{}}}v𝐤|p^\alpha |c𝐤c𝐤|p^\beta |v𝐤`$ (109)
$`\times \delta (E_{c𝐤}E_{v𝐤}\mathrm{}\omega )\text{d}^3k,`$
where $`\mathrm{}\omega `$ is the photon energy, and $`|n𝐤`$ is an eigenket of $`H`$ with energy $`E_{n𝐤}`$. The sum covered the four valence bands $`v`$ and the seven lowest conduction bands $`c`$. The integral was performed using a modified Gilat-Raubenheimer technique Gilat and Raubenheimer (1966); Saravia and Brust (1968); Ren and Harrison (1981) based on 45961 k points in the irreducible part of the Brillouin zone Brust (1964); Cornwell (1984) (representing 2048000 points in the full Brillouin zone). The energy interval for this calculation was 1 meV.

To reveal more clearly the physical meaning of the calculated spectra, the same method was used to calculate the joint density of states function

$$J(E)=\underset{c,v}{}\frac{2\mathrm{\Omega }_0}{(2\pi )^3}_{\mathrm{\Omega }_0^{}}\delta (E_{c𝐤}E_{v𝐤}E)\text{d}^3k.$$
(110)
The dielectric function differs from $`J(E)`$ in that each transition is weighted by the oscillator strength

$$f_{cv}^{\alpha \beta }=\frac{2v𝐤|p^\alpha |c𝐤c𝐤|p^\beta |v𝐤}{m(E_{c𝐤}E_{v𝐤})}.$$
(111)
One may therefore use $`ϵ_2`$ and $`J`$ to define an average oscillator strength at each energy by

$$F^{\alpha \beta }(\mathrm{}\omega )=\frac{m\omega \mathrm{\Omega }_0}{2\mathrm{}\pi ^2e^2}\frac{ϵ_2^{\alpha \beta }(\omega )}{J(\mathrm{}\omega )}.$$
(112)
In cubic crystals, the tensors (109) and (112) reduce to scalars: $`ϵ_2^{\alpha \beta }(\omega )=ϵ_2(\omega )\delta _{\alpha \beta }`$.

The calculated dielectric function $`ϵ_2(\omega )`$ for Ge and Si is plotted in Fig. 3.

This figure compares experimental data Aspnes and Studna (1983) for the dielectric function with the values given by Eq. (109) for (i) the nonlocal pseudopotential model of Chelikowsky and Cohen, Chelikowsky and Cohen (1973, 1974, 1976) and (ii) the 15-orbital tight-binding model of Table 2. For each model, two plots of $`ϵ_2(\omega )`$ are given, corresponding to two different expressions for the momentum operator $`𝐩`$.

In pseudopotential calculations, optical properties are usually calculated from $`𝐀𝐩`$ coupling with $`𝐩=i\mathrm{}`$. Chelikowsky and Cohen (1976) However, if the pseudopotential is nonlocal, this coupling is not gauge invariant; the correct linear coupling is given instead by the kinematic momentum $`𝐩=(m/i\mathrm{})[𝐱,H]`$. Ismail-Beigi et al. (2001) Since pseudopotential calculations are usually performed in a plane-wave basis, a more convenient expression for the kinematic momentum is given by Eq. (88) (which is valid for both discrete and continuous coordinates $`𝐱`$).

In the present tight-binding theory, the kinematic momentum is given by Eq. (108). The two tight-binding functions plotted in Fig. 3 correspond to two different choices of the intra-atomic coordinates $`𝐱_\alpha `$, which are determined by the parameters $`r`$ and $`r^{}`$ in Table 2. One choice was the limit $`r0`$, $`r^{}0`$, which is equivalent to the zero-parameter model of Eqs. (3) and (4). The other, more physically realistic choice was $`r=\frac{1}{3}a`$ and $`r^{}=\sqrt{2}r`$. The value $`r=\frac{1}{3}a`$ was chosen because it yields equidistant lattice sites along the bond directions $`111`$. The value $`r^{}=\sqrt{2}r`$ was used because a somewhat larger value (e.g., $`1.5r`$) breaks the link $`\alpha _{ff}`$ in Table 2, whereas a somewhat smaller value (e.g., $`1.3r`$) breaks the links $`\beta _{bf}`$ and $`\beta _{ff}`$ (while simultaneously forming a new link $`\beta _{be}`$). These values of $`r`$ and $`r^{}`$ also generated successful starting values for some of the parameters in Table 2 (although the final fitted parameters were not very close to the starting values).

Several conclusions may be drawn from Fig. 3. The first is that, within a given model (pseudopotential or tight binding), the choice of momentum operator does not have much numerical significance for the present calculation. This was to be expected on physical grounds, since the intra-atomic coordinate $`𝐱_\alpha `$ (in the tight-binding model) and the nonlocal part of the momentum (in the pseudopotential model) both lead to polarization effects within the atom. Yet it is well known that the bonds between atoms are much easier to polarize than the atoms themselves. Harrison (1989, 1999) Hence, in a bulk semiconductor, intra-atomic effects yield only a minor numerical correction. This conclusion should remain valid in any system where the states are extended, but it may break down in systems where localized states are important. Cruz et al. (1999); Pedersen et al. (2001)

The nonlocal part of the pseudopotential momentum tends to increase $`ϵ_2(\omega )`$ in most frequency ranges, but it sometimes has the opposite effect (see, e.g., the region below 4 eV in Si). However, the intra-atomic coupling in the tight-binding model always decreases $`ϵ_2(\omega )`$. This may be understood by noting that the dominant nonlocal term in the Hamiltonian of Table 2 is the coupling $`\beta _{bb}`$ along the bond between nearest neighbors. Increasing the value of $`r`$ decreases the distance between $`b`$ sites on neighboring atoms, thereby decreasing the momentum matrix in Eq. (108). This tends to increase (slightly) the discrepancy between the tight-binding and pseudopotential dielectric functions. It is possible, however, that a different parametrization of the Hamiltonian might yield different results.

The tight-binding and pseudopotential dielectric functions are quite similar in Ge, but there is a significant discrepancy at the $`E_2`$ peak in Si. The reason for the difference between the models is apparent from the joint density of states $`J`$ and average oscillator strength $`F`$ plotted in Fig. 4.

This figure shows that the tight-binding $`J`$ is very accurate in Ge and somewhat less so in Si, as might have been expected from the quality of the fitted energy bands in Fig. 2. However, the tight-binding model underestimates the oscillator strength over almost the entire frequency range shown, typically by about 20%. The difference is most pronounced between 4 and 4.5 eV in Si, where it exceeds 30%. When combined with a slight underestimate of $`J`$ in the same region, this leads to the discrepancy in the $`E_2`$ peak noted above.

The $`E_2`$ peak in Si is associated with a volume in $`𝐤`$ space near $`(0.9,0.1,0.1)2\pi /a_0`$, Chelikowsky and Cohen (1976) which is close to the $`X`$ point. Thus, the physical reason for the error in the $`E_2`$ peak is probably the spurious $`X_1`$ conduction band near 9 eV in Fig. 2. Because this band is too low in energy, it mixes more strongly with the lowest $`X_1`$ conduction band in the tight-binding model than it does in the pseudopotential model. This change in the wave function causes a corresponding change in oscillator strength. Thus, the majority of the error in the Si $`E_2`$ peak would likely be eliminated if one could find an improved parameter set that raises the energy of the upper $`X_1`$ conduction band.

It is less clear whether the systematic underestimate of oscillator strength at all frequencies could be resolved by changing the Hamiltonian parameters. Oscillator strength was not included in the present fitting routine, so it is possible that with specific attention to this feature, one could improve the oscillator strength while maintaining the quality of the joint density of states. However, it is also possible that such an underestimate is a fundamental limitation imposed by the small basis size in the tight-binding model. Thus, at present, the 15-orbital model is capable of providing semiquantitative predictions of oscillator strength that reproduce all of the major trends exhibited by the pseudopotential model. Whether future developments bring it into precise quantitative agreement remains to be seen.

Finally, it is worth noting that in Fig. 3, the calculated $`E_1`$ peak (between 2 and 3 eV) for Ge is considerably less than the experimental value for both the pseudopotential and tight-binding models. Chelikowsky and Cohen have attributed this discrepancy to the neglect of exciton effects. Chelikowsky and Cohen (1973) However, a recent first-principles calculation of $`ϵ_2(\omega )`$ that includes electron-hole interactions Benedict et al. (1998) shows that the contribution from excitons in Ge is not large enough to fill the gap in Fig. 3. Thus, the empirical pseudopotential for Ge Chelikowsky and Cohen (1973, 1976) probably needs further adjustment to increase the $`E_1`$ peak.

## VI Discussion and conclusions

This paper has shown that intra-atomic optical transitions can be incorporated into tight-binding theory in a gauge-invariant way if the coordinate representation is taken as fundamental. Orthogonal atomic-like orbitals can be constructed from symmetrized coordinate eigenkets, and the coupling to electromagnetic fields can then be described using lattice gauge theory. A model based on 15 such orbitals per atom is capable of describing the most important features of the tetrahedral semiconductors Ge and Si. This basis is slightly larger than existing 10-orbital models, Jancu et al. (1998) but it has the advantages of (i) gauge invariance and (ii) providing an explicit wave function for the electron. A larger basis is needed in the present theory because the restrictions imposed by local gauge symmetry reduce the number of available Hamiltonian fitting parameters.

The field-particle coupling derived here is similar to that given by the Peierls substitution, Peierls (1933); Graf and Vogl (1995); Dumitrică et al. (1998); Boykin et al. (2001); Ismail-Beigi et al. (2001); Boykin and Vogl (2002) but the field-induced phase factor appears in the coordinate representation rather than the tight-binding representation. Thus, the present formalism includes intra-atomic coupling not present in the Peierls substitution. Ismail-Beigi, Chang, and Louie (2001) have recently presented a derivation of the Peierls phase for nonlocal Hamiltonians in a continuous coordinate representation. They have argued that this derivation justifies the use of the Peierls substitution in tight-binding theory. However, their derivation cannot be extrapolated to tight-binding theory, because ordinary $`𝐀𝐩`$ coupling gives rise to intra-atomic interactions that are not included in the Peierls substitution. The existence of such interactions was not considered in the tight-binding theory of Ref. Ismail-Beigi et al., 2001.

It is interesting to consider whether there are any other ways of incorporating local gauge symmetry into tight-binding theory. One possibility is to work directly in the usual tight-binding representation \[see, e.g., Eq. (4)\], where the basis kets are labeled by the symmetry of the orbital and the position of the atom. If gauge symmetry is to be applied in this basis, the coordinate operator must be diagonal, hence all intra-atomic matrix elements must be set to zero \[as in Eq. (3)\]. One could then introduce an Abelian $`U(1)`$ gauge field on this lattice using the approach described above in Sec. IV. The results would be identical to those found in Sec. IV, except that the phase factor in the Hamiltonian (39) would be applied in the tight-binding representation rather than the coordinate representation. Hence, this approach would constitute a “derivation” of the Peierls phase (4). Such a derivation would eliminate the ambiguity associated with the choice of path Graf and Vogl (1995) in Eq. (4). (Other techniques for eliminating path ambiguity are described in Refs. Ismail-Beigi et al., 2001 and Boykin et al., 2001.)

The problem with this approach lies in its treatment of the coordinate operator. In the theory described in Sec. IV, when the basis size is increased, the eigenvalue spectrum of the coordinate operator remains nondegenerate, tending (in the limit of infinite basis dimensions) toward a continuous spectrum. However, if the coordinate operator is required to be diagonal in the tight-binding basis, its eigenvalue spectrum is always degenerate, tending (in the limit of infinite basis dimensions) toward a discrete spectrum with infinite degeneracy. Hence, any tight-binding theory that is either based on or equivalent to the Peierls substitution cannot reproduce the correct continuum limit of the coordinate operator.

As a generalization of the above approach, one might also consider introducing a non-Abelian gauge field Yang and Mills (1954); Wu and Yang (1975); Moriyasu (1983); Guidry (1991); Wilson (1974); Kogut and Susskind (1975); Rebbi (1983); Rothe (1997); Christ et al. (1982a, b, c) in the tight-binding basis. The idea would be to treat the tight-binding electron as a new type of “elementary particle” with some internal degrees of freedom (corresponding to the symmetry labels of the atomic orbitals) that are coupled to the gauge field. In this way, one might hope to reproduce the effects of intra-atomic coupling while remaining in the tight-binding basis. There are, however, numerous difficulties with this approach.

First, the lattice sites in lattice gauge theory represent states of the same particle at different positions. Hence, these states are identical apart from their positions. However, in tight-binding theory the atoms are generally not the same. Second, the field equations for a non-Abelian gauge field are intrinsically nonlinear, because the field carries its own charge and is coupled directly to itself (i.e., it is self-radiating). It is therefore difficult to imagine how such a field could reproduce ordinary electromagnetism, in which the field has no charge, and nonlinearities arise only from interactions with matter. Third, one would need to define a new gauge field theory every time one added new orbitals to the model, and every one of these non-Abelian field theories would need to reproduce the results of Abelian electromagnetic theory. Finally, the coordinate operator in this approach would still have a discrete, degenerate eigenvalue spectrum.

Thus, it appears that the present approach—that is, an Abelian $`U(1)`$ gauge field in the coordinate representation—is the only gauge-invariant method for including electromagnetic fields in empirical tight-binding theory that tends toward the correct continuum limit as the basis dimensions are increased. In this case, the only way that the essential structure of the theory can be modified is to change the topology of the system so as to increase the number of links between lattice sites. This would increase the number of free parameters in the Hamiltonian, thereby permitting a reduction in basis size. Such a modification would clearly be beneficial, but it is not obvious that there exists any alternative topology for general lattices that is capable of reproducing continuum electromagnetism unambiguously. Hence, this possibility will not be explored further here.

###### Acknowledgements.

I am grateful to Peter Vogl, Tim Boykin, and Tai Kai Ng for helpful discussions. This work was supported by Hong Kong RGC Grant No. HKUST6139/00P.

## Appendix A Symmetrized orbitals

This section presents symmetrized orbitals obtained by applying the symmetry operations of the cubic group $`O_h`$ to the coordinate eigenkets $`|100`$ and $`|110`$. \[The orbitals for $`|111`$ may be found in Eq. (II).\] For $`|100`$, if the basis kets are ordered as $`\{|100`$, $`|010`$, $`|001`$, $`|\overline{1}00`$, $`|0\overline{1}0`$, $`|00\overline{1}\}`$, then the symmetrized orbitals are

$`|\mathrm{\Gamma }_1`$ $`=`$ $`|s`$
$`=`$ $`\frac{1}{\sqrt{6}}(1,1,1,1,1,1),`$
$`|\mathrm{\Gamma }_{15}^z`$ $`=`$ $`|p_z`$
$`=`$ $`\frac{1}{\sqrt{2}}(0,0,1,0,0,1),`$
$`|\mathrm{\Gamma }_{12}^a`$ $`=`$ $`|d_{2z^2x^2y^2}`$
$`=`$ $`\frac{1}{2\sqrt{3}}(1,1,2,1,1,2),`$
$`|\mathrm{\Gamma }_{12}^b`$ $`=`$ $`|d_{x^2y^2}`$ (113)
$`=`$ $`\frac{1}{2}(1,1,0,1,1,0).`$
For $`|110`$, if the basis kets are ordered as $`\{|011`$, $`|01\overline{1}`$, $`|0\overline{1}1`$, $`|0\overline{1}\overline{1}`$, $`|101`$, $`|10\overline{1}`$, $`|\overline{1}01`$, $`|\overline{1}0\overline{1}`$, $`|110`$, $`|1\overline{1}0`$, $`|\overline{1}10`$, $`|\overline{1}\overline{1}0\}`$, then the symmetrized orbitals are

$`|\mathrm{\Gamma }_1`$ $`=`$ $`|s`$
$`=`$ $`\frac{1}{2\sqrt{3}}(1,1,1,1,1,1,1,1,1,1,1,1),`$
$`|\mathrm{\Gamma }_{15}^z`$ $`=`$ $`|p_z`$
$`=`$ $`\frac{1}{2\sqrt{2}}(1,1,1,1,1,1,1,1,0,0,0,0),`$
$`|\mathrm{\Gamma }_{12}^a`$ $`=`$ $`|d_{2z^2x^2y^2}`$
$`=`$ $`\frac{1}{2\sqrt{6}}(1,1,1,1,1,1,1,1,2,2,2,2),`$
$`|\mathrm{\Gamma }_{12}^b`$ $`=`$ $`|d_{x^2y^2}`$
$`=`$ $`\frac{1}{2\sqrt{2}}(1,1,1,1,1,1,1,1,0,0,0,0),`$
$`|\mathrm{\Gamma }_{25^{}}^{xy}`$ $`=`$ $`|d_{xy}`$
$`=`$ $`\frac{1}{2}(0,0,0,0,0,0,0,0,1,1,1,1),`$
$`|\mathrm{\Gamma }_{25}^c`$ $`=`$ $`|f_{z(x^2y^2)}`$ (114)
$`=`$ $`\frac{1}{2\sqrt{2}}(1,1,1,1,1,1,1,1,0,0,0,0).`$
For the triply degenerate representations $`\mathrm{\Gamma }_{15}`$, $`\mathrm{\Gamma }_{25^{}}`$, and $`\mathrm{\Gamma }_{25}`$, only one representative orbital is given; the others may be obtained from cyclic permutations of $`x`$, $`y`$, and $`z`$.

## Appendix B Geometry of Voronoi polyhedra

This appendix presents an algorithm for calculating the geometry of the Voronoi polyhedra associated with a given set of nodes $`𝐱_i`$. The basic element in this algorithm is a procedure for finding the edges of the polyhedra. An edge of a Voronoi polyhedron is a finite line segment consisting of points that are closer to three (or more) nodes than to any other nodes. The first step is therefore to determine the equation defining this line.

Any three noncollinear points $`𝐱_i`$, $`𝐱_j`$, and $`𝐱_k`$ define a plane whose normal is the vector

$`𝐧_{ijk}`$ $`=`$ $`𝐝_{ji}\times 𝐝_{ki}`$ (115)
$`=`$ $`𝐱_i\times 𝐱_j+𝐱_j\times 𝐱_k+𝐱_k\times 𝐱_i,`$
where $`𝐝_{ji}=𝐱_j𝐱_i`$. This plane is the set of points $`𝐱`$ satisfying

$$𝐧_{ijk}(𝐱𝐱_i)=0.$$
(116)
The line consisting of all points $`𝐱`$ equidistant from $`𝐱_i`$, $`𝐱_j`$, and $`𝐱_k`$ may therefore be written as

$$𝐱=𝐱_{ijk}+\lambda \widehat{𝐧}_{ijk},$$
(117)
where $`\lambda `$ is a real parameter, $`\widehat{𝐧}_{ijk}=𝐧_{ijk}/n_{ijk}`$, and $`𝐱_{ijk}`$ is the point in the plane (116) equidistant from $`𝐱_i`$, $`𝐱_j`$, and $`𝐱_k`$. To determine this point, note that points $`𝐱`$ equidistant from $`𝐱_i`$ and $`𝐱_j`$ satisfy

$$(𝐱_j𝐱_i)[𝐱\frac{1}{2}(𝐱_j+𝐱_i)]=0.$$
(118)
The point $`𝐱_{ijk}`$ therefore satisfies the three equations

$$𝐚_r𝐱_{ijk}=c_r(r=1,2,3)$$
(119)
in which

$$𝐚_1=𝐝_{ji},𝐚_2=𝐝_{ki},𝐚_3=𝐧_{ijk},$$
(120)
which may be viewed as a set of oblique (more specifically, monoclinic) basis vectors, and

$`c_1`$ $`=`$ $`\frac{1}{2}(x_j^2x_i^2),`$
$`c_2`$ $`=`$ $`\frac{1}{2}(x_k^2x_i^2),`$
$`c_3`$ $`=`$ $`𝐧_{ijk}𝐱_i`$ (121)
$`=`$ $`𝐱_i(𝐱_j\times 𝐱_k).`$
The solution to Eqs. (119) is given by

$$𝐱_{ijk}=\underset{s=1}{\overset{3}{}}c_s𝐛_s,$$
(122)
where $`𝐛_s`$ is a reciprocal basis vector satisfying $`𝐚_r𝐛_s=\delta _{rs}`$; e.g.,

$$𝐛_1=\frac{𝐚_2\times 𝐚_3}{𝐚_1(𝐚_2\times 𝐚_3)}.$$
(123)
Now since $`𝐚_1(𝐚_2\times 𝐚_3)=n_{ijk}^2`$, the vectors $`𝐛_s`$ are given explicitly by

$`𝐛_1`$ $`=`$ $`[𝐝_{ji}(d_{ki}^2)𝐝_{ki}(𝐝_{ji}𝐝_{ki})]/n_{ijk}^2,`$
$`𝐛_2`$ $`=`$ $`[𝐝_{ki}(d_{ji}^2)𝐝_{ji}(𝐝_{ji}𝐝_{ki})]/n_{ijk}^2,`$ (124)
$`𝐛_3`$ $`=`$ $`𝐧_{ijk}/n_{ijk}^2.`$
Equation (122) may then be rearranged in the more symmetric form

$$𝐱_{ijk}=\frac{𝐱_i[d_{jk}^2(𝐝_{ji}𝐝_{ki})]+𝐱_j[d_{ki}^2(𝐝_{kj}𝐝_{ij})]+𝐱_k[d_{ij}^2(𝐝_{ik}𝐝_{jk})]}{2n_{ijk}^2}.$$
(125)
This result, together with Eq. (117), defines the line equidistant from nodes $`𝐱_i`$, $`𝐱_j`$, and $`𝐱_k`$.

The next step is to determine whether any segment of this line forms an edge of a Voronoi polyhedron. Points on such a segment must lie closer to $`𝐱_i`$, $`𝐱_j`$, and $`𝐱_k`$ than to any other node $`𝐱_l`$. For each node $`𝐱_l`$, one calculates

$$\alpha _l=𝐝_{il}\widehat{𝐧}_{ijk}.$$
(126)
If $`\alpha _l=0`$, then $`𝐱_l`$ lies in the plane (116). In this case, if $`|𝐱_l𝐱_{ijk}|<|𝐱_i𝐱_{ijk}|`$, then no portion of the line (117) forms an edge of a Voronoi polyhedron. On the other hand, if $`|𝐱_l𝐱_{ijk}||𝐱_i𝐱_{ijk}|`$, then the line (117) may form an edge (depending on the position of the other nodes $`𝐱_l^{}`$).

If $`\alpha _l0`$, then $`𝐱_l`$ does not lie in the plane (116). In this case, points on the line (117) that are closer to $`𝐱_i`$ than to $`𝐱_l`$ satisfy \[cf. Eq. (118)\]

$$𝐝_{il}(𝐱_{ijk}+\lambda \widehat{𝐧}_{ijk}𝐱_{il})>0,$$
(127)
in which $`𝐱_{il}=\frac{1}{2}(𝐱_i+𝐱_l)`$. Hence, the position of the point equidistant from $`𝐱_i`$, $`𝐱_j`$, $`𝐱_k`$, and $`𝐱_l`$ is given by the following value of the parameter $`\lambda `$ in Eq. (117):

$$\lambda _l=\frac{𝐝_{il}(𝐱_{il}𝐱_{ijk})}{\alpha _l}.$$
(128)
One may then define

$$\lambda _{\text{min}}=\text{max}(\lambda _l|\alpha _l>0)$$
(129)
(i.e., the maximum value of $`\lambda _l`$ for all $`l`$ such that $`\alpha _l>0`$) and

$$\lambda _{\text{max}}=\text{min}(\lambda _l|\alpha _l<0).$$
(130)
If $`\lambda _{\text{max}}>\lambda _{\text{min}}`$, then the line segment (117) with $`\lambda _{\text{min}}<\lambda <\lambda _{\text{max}}`$ forms an edge of a Voronoi polyhedron. This establishes the positions of two corners of the polyhedron:

$`𝐱_c`$ $`=`$ $`𝐱_{ijk}+\lambda _{\text{min}}\widehat{𝐧}_{ijk},`$
$`𝐱_c^{}`$ $`=`$ $`𝐱_{ijk}+\lambda _{\text{max}}\widehat{𝐧}_{ijk}.`$ (131)
The set of all nodes in the plane (116) that lie closer to the line segment $`\lambda _{\text{min}}<\lambda <\lambda _{\text{max}}`$ than any other node defines what is called a plaquette. Since there are in general more than three such nodes, it is convenient to define a unique label $`q`$ for each plaquette, with

$$\widehat{𝐧}_q\widehat{𝐧}_{ijk},𝐱_q𝐱_{ijk},$$
(132)
for any members $`𝐱_i`$, $`𝐱_j`$, and $`𝐱_k`$ of the given plaquette. (The sign of $`\widehat{𝐧}_q`$ is fixed by some convention for the ordering of the nodes $`𝐱_i`$, $`𝐱_j`$, and $`𝐱_k`$.) Each plaquette is associated uniquely with one edge of a Voronoi polyhedron, the length of which is

$$d_q=|𝐱_c^{}𝐱_c|=\lambda _{\text{max}}\lambda _{\text{min}},$$
(133)
with $`d_q>0`$ by definition.

At this point, one has sufficient information to determine whether a link exists between any pair of nodes $`𝐱_i`$ and $`𝐱_j`$. The first step is to use the above procedure to find all of the corner points $`𝐱_c`$ common to nodes $`i`$ and $`j`$. By definition, all such points lie in the plane (118). The set of these points defines a polygon, the perimeter of which consists of the line segments (133). The area of the polygon may be calculated by numbering the corner points $`𝐱_c`$ in sequential order around the perimeter of the polygon, then partitioning the polygon into triangles as shown in Fig. 5.

The normal area vector is

$$𝐒_{ij}=\frac{1}{2}\underset{c=3}{\overset{N_{ij}}{}}(𝐱_{c1}𝐱_1)\times (𝐱_c𝐱_1),$$
(134)
where $`N_{ij}`$ is the number of corner points common to nodes $`i`$ and $`j`$. The area $`S_{ij}=|𝐒_{ij}|`$ is the area of the surface shared by the Voronoi polyhedra for sites $`i`$ and $`j`$; note that if $`N_{ij}=1`$ or 2, the shared region is a point or line, and the area (134) is zero. Nodes $`𝐱_i`$ and $`𝐱_j`$ are linked only if $`S_{ij}>0`$.

The volume $`\mathrm{\Omega }_i`$ of the Voronoi polyhedron for node $`𝐱_i`$ may be calculated from $`S_{ij}`$ and $`d_{ij}`$. One simply integrates the identity $`𝐱=3`$ over the polyhedron, using the divergence theorem and the fact that the plane containing $`S_{ij}`$ is the perpendicular bisector of $`𝐝_{ij}`$ (although $`𝐝_{ij}`$ need not intersect $`S_{ij}`$ itself). The result is

$$\mathrm{\Omega }_i=\frac{1}{6}\underset{j}{}S_{ij}d_{ij}.$$
(135)
For each link $`\mathrm{}(i,j)`$, one can construct a polyhedron by drawing lines from nodes $`𝐱_i`$ and $`𝐱_j`$ to each of their common corner points $`𝐱_c`$. The volume of this polyhedron is

$$\mathrm{\Omega }_{\mathrm{}}=\frac{1}{3}S_{\mathrm{}}d_{\mathrm{}}.$$
(136)
The volume $`\mathrm{\Omega }_{\mathrm{}}\mathrm{\Omega }_{ij}`$ is bisected by $`S_{\mathrm{}}S_{ij}`$, with half lying in $`\mathrm{\Omega }_i`$ and half in $`\mathrm{\Omega }_j`$; hence

$$\mathrm{\Omega }_i=\frac{1}{2}\underset{j}{}\mathrm{\Omega }_{ij}.$$
(137)
The nodes in plaquette $`q`$ define a polygon in the plane (116); the perimeter of this polygon is formed by the links $`d_{\mathrm{}}`$. Hence, the area of the plaquette can be calculated in the same way as the link area (134):

$$𝐒_q=\frac{1}{2}\underset{i=3}{\overset{N_q}{}}(𝐱_{i1}𝐱_1)\times (𝐱_i𝐱_1),$$
(138)
where $`N_q`$ is the number of nodes in plaquette $`q`$. A polyhedron may be constructed for each plaquette by drawing lines from each node of the plaquette to the corner points (131); the volume of this polyhedron is

$$\mathrm{\Omega }_q=\frac{1}{3}S_qd_q.$$
(139)
Finally, the plaquette surfaces $`S_q`$ partition all of space into nonoverlapping polyhedra (this is referred to as a Delaunay tessellation Okabe et al. (2000)). These polyhedra (or cells) are in one-to-one correspondence with the corner points $`𝐱_c`$ of the Voronoi polyhedra. The volume of cell $`c`$ is

$$\mathrm{\Omega }_c=\frac{1}{3}\underset{qc}{}𝐒_q(𝐱_q𝐱_c),$$
(140)
where the direction of $`𝐒_q`$ is chosen to point outward from $`\mathrm{\Omega }_c`$ (note that $`𝐱_c`$ does not necessarily lie inside $`\mathrm{\Omega }_c`$, Christ et al. (1982a) so the dot product may be negative for some $`q`$).

A useful set of sum rules for verifying the consistency of a calculated geometry is

$$\underset{i}{}\mathrm{\Omega }_i=\underset{\mathrm{}}{}\mathrm{\Omega }_{\mathrm{}}=\underset{q}{}\mathrm{\Omega }_q=\underset{c}{}\mathrm{\Omega }_c=\mathrm{\Omega },$$
(141)
where $`\mathrm{\Omega }`$ is the volume of some region over which the node distribution is periodic, such as a primitive cell in a Bravais lattice (not to be confused with the generally nonperiodic cell $`\mathrm{\Omega }_c`$). The sum rule for $`\mathrm{\Omega }_i`$ follows directly from the definition of $`\mathrm{\Omega }_i`$ given in Sec. III, since every point in $`\mathrm{\Omega }`$ must lie in at least one Voronoi polyhedron, and the only regions of overlap between polyhedra are points, lines, or planes of zero volume. The sum rule for $`\mathrm{\Omega }_c`$ was proven in Ref. Christ et al., 1982a. The sum rule for $`\mathrm{\Omega }_{\mathrm{}}`$ follows from that for $`\mathrm{\Omega }_i`$, since the set $`\{\mathrm{\Omega }_{\mathrm{}}\}`$ is just another way of partitioning the set $`\{\mathrm{\Omega }_i\}`$ \[see Eq. (137)\]. Likewise, the sum rule for $`\mathrm{\Omega }_q`$ follows from that for $`\mathrm{\Omega }_c`$, since the set $`\{\mathrm{\Omega }_q\}`$ is just another way of partitioning the set $`\{\mathrm{\Omega }_c\}`$ \[see Eqs. (133), (139), and (140)\].

## Appendix C Examples of link geometry

This appendix presents values of the link lengths $`d_{\mathrm{}}`$ and surface areas $`S_{\mathrm{}}`$ for several lattices. The simplest geometry occurs for Bravais lattices, of which only the cubic lattices are considered here. For the simple cubic lattice, only nearest neighbors are linked, with $`d_1=a`$ and $`S_1=a^2`$, where $`a`$ is the lattice constant. For the body-centered cubic lattice, both first and second nearest neighbors are linked, with

$$d_1=\frac{\sqrt{3}}{2}a,S_1=\frac{3\sqrt{3}}{16}a^2,$$
(142)
and

$$d_2=a,S_2=\frac{1}{8}a^2.$$
(143)
For the face-centered cubic lattice, only nearest neighbors are linked, with

$$d_1=\frac{a}{\sqrt{2}},S_1=\frac{a^2}{4\sqrt{2}}.$$
(144)
The remaining lattices to be considered are those obtained by putting symmetrized orbitals on the atomic sites of the diamond or zinc-blende structure. If only a single $`s`$ orbital per atom is used (i.e., one $`|000`$ basis ket per atom), then each atom is linked to 4 nearest neighbors and 12 second-nearest neighbors, with

$$d_1=\sqrt{3}a,S_1=3\sqrt{3}a^2,$$
(145)
and

$$d_2=2\sqrt{2}a,S_2=\frac{\sqrt{2}}{4}a^2.$$
(146)
Here $`a=\frac{1}{4}a_0`$, where $`a_0`$ is the conventional cubic lattice constant. Note that in this case, the link $`d_2`$ does not intersect the surface $`S_2`$.

Coupling between second-nearest neighbors persists in models with more than one orbital per atom. In the “$`sp^3`$” model with four $`|111`$ sites per atom (generated by applying the symmetry operations of $`T_d`$ to $`|r,r,r`$ and $`|ar,ar,ar`$), each site is linked to three others on the same atom:

$$d_0=2\sqrt{2}r,S_0=\frac{7\sqrt{2}}{4}a^2,$$
(147)
one on a neighboring atom:

$$d_1=\sqrt{3}(a2r),S_1=3\sqrt{3}a^2,$$
(148)
and three second-nearest neighbors:

$$d_2=2\sqrt{2}(ar),S_2=\frac{\sqrt{2}}{4}a^2.$$
(149)
In the model generated by either $`T_d`$ or $`O_h`$ and $`|r,0,0`$, each of the six sites is linked to four others on the same atom:

$$d_0=\sqrt{2}r,S_0=\frac{a^2(5a6r)}{4\sqrt{2}(ar)},$$
(150)
four nearest neighbors:

$$d_1=\sqrt{3a^24ar+2r^2},S_1=\frac{a^2d_1}{2(ar)},$$
(151)
and four second-nearest neighbors:

$$d_2=\sqrt{2}(2ar),S_2=\frac{a^2(a2r)}{4\sqrt{2}(ar)}.$$
(152)
In the model generated by either $`T_d`$ or $`O_h`$ and $`|r,r,0`$, each of the 12 sites is linked to four others on the same atom, two of which have

$$d_0=\sqrt{2}r,S_0^{}=\frac{a^2(3a4r)}{4\sqrt{2}(ar)},$$
(153)
and two of which have

$$d_0=\sqrt{2}r,S_0^{\prime \prime }=\frac{a^2(7a8r)}{4\sqrt{2}(ar)}.$$
(154)
Each site is also linked to two nearest neighbors:

$$d_1=\sqrt{3a^28ar+6r^2},S_1=\frac{a^2d_1}{2(ar)},$$
(155)
and one second-nearest neighbor:

$$d_2=2\sqrt{2}(ar),S_2=\frac{a^3}{2\sqrt{2}(ar)}.$$
(156)
In the model generated by $`O_h`$ and $`|r,r,r`$, there are two distinct lattice sites. Sites such as $`|r,r,r`$ and $`|ar,ar,ar`$ are labeled $`b`$ because they lie on the bonds between atoms, whereas sites such as $`|r,r,r`$ and $`|a+r,a+r,a+r`$ are labeled $`e`$ because they lie on “empty” bonds. (Both of these sites are actually Wyckoff $`e`$ sites, but they are inequivalent, because the site symmetry of the atoms in diamond is $`T_d`$.) Each $`b`$ site is linked to three $`e`$ sites on the same atom (and vice versa), with

$$d_0^{be}=2r,S_0^{be}=\frac{a(9a4r)}{8}.$$
(157)
Each $`b`$ site is linked to one nearest-neighbor $`b`$ site:

$$d_1^{bb}=\sqrt{3}(a2r),S_1^{bb}=\frac{3\sqrt{3}}{4}a^2,$$
(158)
and three nearest-neighbor $`e`$ sites:

$$d_1^{be}=\sqrt{a^2+2(a2r)^2},S_1^{be}=\frac{1}{8}ad_1^{be}.$$
(159)
Each $`e`$ site is also linked to three $`b`$ sites via (159), to six nearest-neighbor $`e`$ sites via

$$d_1^{ee}=\sqrt{2a^2+(a2r)^2},S_1^{ee}=\frac{1}{4}ad_1^{ee},$$
(160)
and to three second-nearest neighbor $`e`$ sites via

$$d_2^{ee}=2\sqrt{2}(ar),S_2^{ee}=\frac{\sqrt{2}}{4}a(a+2r).$$
(161)