text:
# First-principles envelope-function theory for lattice-matched semiconductor heterostructures

## I Introduction

### I.1 Background and motivation

Envelope-function models continue to play a key role in the design and interpretation of experiments on semiconductor heterostructures. The canonical “envelope-function approximation,” which by definition makes use of only bulk effective-mass parameters and heterojunction band offsets, has been very successful in explaining a wide range of experiments. Bastard (1988); Bastard et al. (1991) However, recent studies have increasingly emphasized interface-related effects lying outside the scope of conventional envelope-function theory, such as optical Krebs and Voisin (1996); Krebs et al. (1998); Platonov et al. (1999) and electrical Takhtamirov and Volkov (2000a); Reker et al. (2002); Rössler and Kainz (2002) anisotropy, intervalley mixing, Pulsford et al. (1989); Aleĭner and Ivchenko (1993); Fu et al. (1993); Im et al. (1999, 2000); Klipstein (2001); Glinskii et al. (2000) and spin polarization phenomena. Ganichev et al. (2001); Ganichev et al. (2002a, b); Golub (2003); Ganichev et al. (2003); Bel’kov et al. (2003); Ganichev and Prettl (2003); Ganichev et al. (2004); Sinova et al. (2004); Shen et al. (2004); Zhang and Yang (2005); Bernevig and Zhang (2005) These effects are generated (wholly or in part) by interface band-mixing terms in the heterostructure Hamiltonian, such as valence-band mixing of light and heavy holes, Aleĭner and Ivchenko (1992); Ivchenko et al. (1993, 1996) $`\mathrm{\Gamma }`$$`X`$ coupling, Aleĭner and Ivchenko (1993); Fu et al. (1993) and the conduction or valence band Rashba coupling. Winkler (2003)

A great deal of progress in the modeling of these effects can be made on the basis of symmetry information alone. The method of invariants, which was originally developed for bulk semiconductors, Luttinger (1956); Suzuki and Hensel (1974); Bir and Pikus (1974); Trebin et al. (1979) has proved a powerful tool in the study of heterostructures as well. Ivchenko and Pikus (1997); Aleĭner and Ivchenko (1992); Ivchenko et al. (1993, 1996); Glinskii and Kravchenko (1998); Glinskiĭ and Kravchenko (1998); Glinskii et al. (2000); Takhtamirov and Volkov (2000a); Ganichev et al. (2001); Ganichev et al. (2002a, b); Golub (2003); Ganichev et al. (2003); Bel’kov et al. (2003); Ganichev and Prettl (2003); Ganichev et al. (2004); Sinova et al. (2004); Shen et al. (2004); Zhang and Yang (2005); Bernevig and Zhang (2005); Winkler (2003) The standard method of invariants uses symmetry information to construct an explicit interface Hamiltonian, but the same information (supplemented by hermiticity or current-conservation requirements) may also be used to construct connection rules for the envelope functions on opposite sides of the interface. Laikhtman (1992); de-Leon et al. (1998); Kisin (1993, 1994); Kisin et al. (1998); Tokatly et al. (2002); Rodina et al. (2002, 2003) In either case, one obtains a phenomenological model containing some interface parameters whose value is determined by comparison with experiment.

Useful as this approach may be, it does not provide any information about the magnitude of the interface parameters, some of which (for reasons unrelated to symmetry) may happen to be zero or negligibly small. To obtain this information (along with a deeper understanding of the physical origin of the interface phenomena), one must turn to a more detailed microscopic model of the interface. Thus, numerous envelope-function models have been derived directly from the microscopic potential energy; these include interface Hamiltonians Leibler (1975, 1977); Burt (1988, 1992, 1994, 1999); Young (1989); Cuypers and van Haeringen (1991); Karavaev and Tikhodeev (1991); Karavaev and Krivorotov (1996a, b); Foreman (1993, 1995a, 1996, 1998, 2001); Burt and Foreman (1999); Klipstein (1999); Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b); Cortez et al. (2000, 2001) and connection rules Smith and Mailhiot (1986, 1990); Mailhiot and Smith (1990); Trzeciakowski (1988); Ando et al. (1989); Cuypers and van Haeringen (1992); Cuypers and van Haeringen (1993a, b); Aversa and Sipe (1993, 1994); Grinyaev et al. (1994) as well as numerical approaches based on solving the $`𝐤𝐩`$ equations in momentum space. Wang and Zunger (1996); Yi and Razeghi (1997)

However, all of the cited derivations are based on empirical pseudopotentials, in which the potential energy is determined by the choice of some specific model for the interface. This leads to ambiguity in the results, as different choices may yield conflicting predictions as to which interface parameters are important, Foreman (1998); Klipstein (1999); Im et al. (1999, 2000); Klipstein (2001); Takhtamirov and Volkov (2000b) or whether interface band-mixing effects are directly related to band offsets. Foreman (1999); Cortez et al. (2000, 2001); Magri and Zunger (2000) Such conflicts can only be resolved by deriving the envelope-function Hamiltonian from an ab initio self-consistent potential, as suggested by Sham and Lu. Sham and Lu (1989) The purpose of this paper is to present such a derivation and examine its implications for interface band-mixing effects. Numerical applications of the theory are not considered here. A preliminary account of these results has been presented elsewhere. Foreman (2005a)

### I.2 Basic assumptions and limitations

The current “standard model” for condensed-matter physics is based on density-functional theory in the local density approximation (LDA). Hohenberg and Kohn (1964); Kohn and Sham (1965); Payne et al. (1992) However, it is well known that this ground-state formalism does not accurately predict energy gaps in semiconductors. Thus, the present state-of-the-art in band-structure theory involves calculations of the electron self-energy in the shielded-interaction Baym and Kadanoff (1961); Baym (1962) or $`GW`$ approximation. Hedin (1965); Hedin and Lundqvist (1969); Aulbur et al. (2000); Onida et al. (2002) This field is not yet fully mature, as apparent early successes Hybertsen and Louie (1986a); Godby et al. (1988); Zhu and Louie (1991) have been questioned in light of recent developments. Ummels et al. (1998); Schöne and Eguiluz (1998); Arnaud and Alouani (2000); Ku and Eguiluz (2002); Lebègue et al. (2003) Nevertheless, any future refinements in approximation techniques for calculating the quasiparticle band structure will continue to be based on the self-energy in Dyson’s equation. For this reason, the present work relies only upon general properties of the self-energy operator, not upon any specific approximation to this quantity.

For simplicity, however, this paper assumes that the ionic cores are described in terms of norm-conserving pseudopotentials. Hamann et al. (1979); Kleinman (1980); Bachelet and Schlüter (1982); Bachelet et al. (1982); Troullier and Martins (1991) The projector augmented wave (PAW) method of Blöchl, Blöchl (1994); Arnaud and Alouani (2000); Lebègue et al. (2003) which does not involve this approximation, has become increasingly popular in recent years. A generalization of the present theory to incorporate the PAW formalism should be possible, but to avoid undue complexity it is not considered here. In addition, this paper considers only lattice-matched systems (with no bulk or interface strain).

The key assumption used in the present work is that the heterostructure is a small perturbation Andreoni et al. (1978) with respect to some periodic virtual reference crystal (such as Al0.5Ga0.5As for a GaAs/AlAs heterostructure). This assumption has been verified within LDA for both isovalent and heterovalent systems, including GaAs/AlAs, Resta et al. (1989); Baroni et al. (1989) Ge/GaAs, Baroni et al. (1988, 1989) In0.53Ga0.47As/InP, Baroni et al. (1989); Peressi et al. (1990) GaAs/Si/AlAs and GaAs/Ge/AlAs, Peressi et al. (1991); Baroni et al. (1992) Si/Ge, Colombo et al. (1991); Peressi and Baroni (1994) and InAs/GaSb. Peressi et al. (1996); Montanari et al. (1996) Within the framework of pseudopotential theory, one would expect it to be valid for other similar heterostructures also. This assumption makes it possible to describe the self-consistent heterostructure perturbation in terms of nonlinear response theory, with the linear response predominant and the quadratic response providing a weak correction. Resta et al. (1989); Baroni et al. (1988, 1989); Peressi et al. (1990, 1991); Colombo et al. (1991); Baroni et al. (1992); Peressi and Baroni (1994); Peressi et al. (1996); Montanari et al. (1996)

Such an approach has been used by Sham Sham (1966) to derive an effective-mass equation for shallow impurity states in semiconductors. Sham’s work was extended in the preceding paper Foreman (2005b) to obtain expressions for the self-energy at small values of the crystal momentum in lattice-matched heterostructures described by spin-dependent nonlocal pseudopotentials. These expressions are used here to construct a multi-band effective-mass theory for lattice-matched heterostructures.

Although the basic formalism developed here is quite general, to keep the paper to a reasonable length it is necessary to impose some restrictions on the material systems that are treated in detail. (This allows the power series in Sec. V to be terminated at a reasonably low order.) Since only lattice-matched systems are considered here, it is assumed that the symmetry of the reference crystal does not support a spontaneous polarization, because that would generate macroscopic electric fields and piezoelectric strain fields inconsistent with the lattice-matching assumption. Therefore, the present theory is not directly applicable to wurtzite materials.

In addition, it is assumed that the interfaces in a heterovalent system such as Ge/GaAs are macroscopically neutral, Baroni et al. (1989); Harrison et al. (1978) so that the atoms can be grouped together into neutral clusters (of fractional atoms) that carry a dipole moment only at the interface. Foreman (2005b); not (a) Thus, the main application of the present theory would be to isovalent or heterovalent heterostructures of semiconductors with the zinc-blende or diamond structure. However, the only examples treated explicitly here are from isovalent systems.

Lattice mismatch could not be included in this theory without fundamental changes to account for atomic relaxation. Nevertheless, the present results provide a solid foundation for subsequent extensions of the theory to include this effect. Other restrictions, such as macroscopic neutrality and the lack of bulk dipole terms, could be lifted merely by extending some power series expansions to higher order (although, depending on the accuracy that is desired, not (a) this may require the inclusion of certain nonanalytic terms neglected in Ref. Foreman, 2005b).

The fact that (within the pseudopotential approximation) a typical heterostructure is a weak perturbation Andreoni et al. (1978); Resta et al. (1989); Baroni et al. (1988, 1989); Peressi et al. (1990, 1991); Colombo et al. (1991); Baroni et al. (1992); Peressi and Baroni (1994); Peressi et al. (1996); Montanari et al. (1996) makes possible the existence of energy eigenstates $`|\psi `$ whose wave function $`n𝐤|\psi `$ in the Luttinger-Kohn representation Luttinger and Kohn (1955) $`|n𝐤`$ of the reference crystal is negligible outside a small region in $`𝐤`$ space near the high-symmetry points of the Brillouin zone. Xia (1989); Wang et al. (1997); Wang and Zunger (1997, 1999) These eigenstates are just the low-energy excitations that are of greatest interest experimentally. Their existence, which has been verified by extensive numerical work on empirical pseudopotential models, Xia (1989); Wang et al. (1997); Wang and Zunger (1997, 1999) is precisely the condition needed for the validity of an effective-mass approximation for slowly varying envelope functions $`F_n(𝐱)`$, Luttinger and Kohn (1955); Leibler (1975, 1977); Burt (1988, 1992, 1994, 1999) since $`F_n(𝐤)n𝐤|\psi `$. The existence of slowly varying envelopes Burt (1992, 1994); Xia (1989); Wang et al. (1997); Wang and Zunger (1997, 1999) provides the foundation for the effective-mass theory derived here.

The present approach, based on quadratic response theory, Foreman (2005b) fits well with the $`𝐤𝐩`$ perturbation formalism of Leibler, Leibler (1975, 1977); Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b) in which the heterostructure perturbation is treated as small in comparison to the energy separating the bands of interest from other remote bands. Leibler’s theory is used here to develop a multi-band effective-mass theory that includes all terms of the same order as the position dependence of the effective mass. This includes cubic and quartic dispersion terms in the Hamiltonian of the reference crystal, as well as the leading contributions from the quadratic response. As shown in Sec. II, higher-order terms cannot consistently be described in terms of local differential equations. Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b)

### I.3 Summary of key results

In this paper, it is shown that the dominant interface band-mixing terms are those arising from the linear response. Indeed, the linear response is the only contribution to band mixing that fits within the framework of the perturbation theory defined above, and the only contribution that can consistently be described in terms of local differential equations. To this level of accuracy, the quadratic response contributes only a bulk band offset term and an interface dipole term, neither of which produces any band mixing in the effective-mass Hamiltonian.

This represents a major simplification, since it implies that the interface band-mixing terms in the Hamiltonian are just a superposition of parameters derived from the linear response to individual ionic perturbations (or neutral cluster perturbations Foreman (2005b); not (a) in the case of heterovalent substitutions). These parameters are calculated once and for all for a given material system; they then appear as coefficients in front of structure factors describing the distribution of atoms in the heterostructure. The Hamiltonian contains spatial derivatives of (and differences between) these structure factors, which generate $`\delta `$-like functions at an interface. The atomic distribution functions change for heterostructures of different shape (wells, wires, or dots) and orientation, but their coefficients do not. Thus, a single set of interface parameters governs the band mixing for any type of heterostructure not (b) in a given material system.

If desired, one can express these linear parameters as a difference between the properties of the various bulk materials that make up the heterostructure. (A similar result was obtained for linear band offsets in Refs. Resta et al., 1989; Baroni et al., 1988, 1989; Peressi et al., 1990, 1991; Colombo et al., 1991; Baroni et al., 1992; Peressi and Baroni, 1994; Peressi et al., 1996; Montanari et al., 1996.) However, in a no-common-atom system such as InAs/GaSb, this must include the (lattice-matched) “interface” materials GaAs and InSb. Also, these bulk-like linear response parameters (which in principle do not require a supercell calculation Baroni et al. (2001)) cannot be determined from experiments on bulk semiconductors.

For some terms in the linear response (derived from the analytic part of the self-energyForeman (2005b)), the interface Hamiltonian itself has the form of a macroscopic average of the Dirac $`\delta `$ function or its derivatives. For the remaining terms (derived from the nonanalytic part of the self-energyForeman (2005b)), the effective charge density has this localized form, but the potential energy is not as well localized, having the form of a long-range multipole potential. However, in a quantum well or any other heterostructure with two-dimensional translation symmetry, all terms are well localized except for the interface dipole terms (which merely modify the band offsets).

The band-mixing Hamiltonian has the same general form for slowly graded structures (within the virtual crystal approximation) and abrupt heterojunctions, the only difference being the rate of change of the structure factor. That is, in contrast to previous theories based on model potentials or empirical pseudopotentials, no new band-mixing parameters appear at an abrupt junction. Indeed, the localized interface terms in the Hamiltonian derived here are qualitatively identical to the Hamiltonian derived by Leibler Leibler (1975, 1977) for slowly graded heterostructures.

These results shed light on recent suggestions Klipstein (1999); Im et al. (1999, 2000); Klipstein (2001); Takhtamirov and Volkov (2000b) that the $`\mathrm{\Gamma }_1`$$`X_{1z}`$ and $`X_{1x}`$$`X_{1y}`$ intervalley mixing potentials at an ideal (001) GaAs/AlAs heterojunction should be proportional to $`\delta (z)`$. It is shown here that such mixing arises only from the quadratic response and is therefore negligible in comparison to $`\mathrm{\Gamma }_1`$$`X_{3z}`$, $`X_{1z}`$$`X_{3z}`$, and $`X_{3x}`$$`X_{3y}`$ mixing. There is, however, a linear contribution to the $`\mathrm{\Gamma }_1`$$`X_{1z}`$ and $`X_{1x}`$$`X_{1y}`$ mixing that is proportional to $`\delta ^{}(z)`$, as well as a linear $`\delta (z)`$ mixing for nonideal interfaces or when spin-orbit coupling is included. These contributions may help to explain the experimental observations in Refs. Im et al., 1999, 2000; Klipstein, 2001.

### I.4 Outline of the paper

The paper begins in Sec. II with a discussion of which terms are to be retained in $`𝐤𝐩`$ perturbation theory, based upon a review of important recent work by Takhtamirov and Volkov. Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b) The main results of the quadratic response theory developed in the preceding paper Foreman (2005b) are presented (in modified form) in Sec. III. The basic envelope-function formalism is developed in Sec. IV, while Sec. V describes the power series expansions that are used to obtain approximate expressions valid for slowly varying envelopes. Perturbation theory is used to eliminate the coupling to “remote” bands in Sec. VI, yielding the basic expression for the multi-band envelope-function Hamiltonian. Modifications to the Hamiltonian that are necessary if one wishes to describe the material parameters of an abrupt heterojunction as piecewise constant are discussed in Sec. VII. Symmetry properties of the Hamiltonian are discussed in Sec. VIII, where explicit matrix representations of the material properties are given for semiconductors with the zinc-blende structure. Finally, the significance of the results and their relation to previous work in the literature are discussed in Sec. IX.

## II Order of terms included

Takhtamirov and Volkov Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b) have recently demonstrated (using an instructive analogy to the leading relativistic corrections to the Schrödinger equation) that most derivations of heterostructure effective-mass equations in the literature (including those of the present author) do not consistently include all perturbative corrections of the same order. This section presents a review and discussion of their results, with the objective of establishing which terms are to be retained in subsequent perturbative approximations.

The first case to be considered is a wide-gap system satisfying

$$\mathrm{\Delta }\overline{V}\overline{E}_\text{g},$$
(1)
where $`\mathrm{\Delta }\overline{V}`$ is a typical heterojunction band offset, and $`\overline{E}_\text{g}`$ is a typical energy gap (for the virtual bulk reference crystal) separating the band in question from all other “remote” bands. It is assumed that in the reference crystal this band is describable by an effective-mass equation with effective mass $`m^{}`$. This will be the case if the $`𝐤𝐩`$ interaction with remote bands

$$\frac{\mathrm{}\overline{k}\overline{p}}{m}(\overline{k}\lambda )\overline{E}_\text{g}$$
(2)
can be treated as a small perturbation. Here $`m`$ is the free-electron mass, $`\overline{k}2\pi /L`$ is a typical envelope function wavenumber for a quantum well of width $`L`$, $`\overline{p}2\pi \mathrm{}/a`$ is a typical interband momentum matrix element (where $`a`$ is the lattice constant), and it is assumed that $`La`$. Equation (2) defines a length parameter $`\lambda =\mathrm{}\overline{p}/m\overline{E}_\text{g}`$, which (if the free-electron contribution to $`m^{}`$ is negligible) may also be written as $`\lambda \mathrm{}(2m^{}\overline{E}_\text{g})^{1/2}`$. For GaAs, both expressions for $`\lambda `$ give $`\lambda 6`$ Å, thus $`\lambda a`$ and $`\overline{k}\lambda 1`$ in a wide quantum well. For the remainder of this paper the parameters $`\lambda `$ and $`a`$ are used interchangeably, although in general error estimates should be based on the larger of the two quantities.

In general one is interested in cases where the kinetic energy is comparable to the band offset:

$$\frac{\mathrm{}^2\overline{k}^2}{2m^{}}\mathrm{\Delta }\overline{V}(\overline{k}\lambda )^2\overline{E}_\text{g}.$$
(3)
This is the order of terms included in ordinary effective-mass theory, Luttinger and Kohn (1955) which may be compared with the nonrelativistic (Schrödinger) approximation to the Dirac equation. The position dependence of the effective mass can be calculated by treating $`\mathrm{\Delta }\overline{V}/\overline{E}_\text{g}`$ as a perturbation, Leibler (1975, 1977) yielding a correction of order

$$\frac{\mathrm{}^2\overline{k}^2}{2m^{}}\left(\frac{\mathrm{\Delta }\overline{V}}{\overline{E}_\text{g}}\right)\frac{(\mathrm{\Delta }\overline{V})^2}{\overline{E}_\text{g}}(\overline{k}\lambda )^2\mathrm{\Delta }\overline{V},$$
(4)
which is of the same order as the nonparabolic $`k^4`$ terms in the kinetic energy of the reference crystal. Hence, these terms (analogous to the relativistic mass corrections to the Schrödinger equation) must be included if one is to retain all terms of the same order as the position dependence of the effective mass.

As will be seen below, the interface band-mixing terms proportional to a $`\delta `$ function are of order

$$A\delta (z)(\overline{k}a)\mathrm{\Delta }\overline{V},$$
(5)
while those proportional to the derivative of a $`\delta `$ function are of order

$$B\delta ^{}(z)(\overline{k}a)^2\mathrm{\Delta }\overline{V},$$
(6)
in which $`A`$ and $`B`$ are constants. Hence, the latter terms are also comparable to the position dependence of the effective mass.

For a step-function discontinuity $`\mathrm{\Delta }V`$ in the potential energy, the envelope function $`F`$ has a discontinuity in its second derivative given by

$$\frac{\mathrm{\Delta }F^{\prime \prime }}{F}=\frac{2m^{}}{\mathrm{}^2}\mathrm{\Delta }V,$$
(7)
where $`FL^{1/2}`$. This gives rise to an asymptotic behavior in $`k`$ space of $`F(k)(\mathrm{\Delta }F^{\prime \prime })k^3`$. But in the exact Luttinger-Kohn envelope-function representation, Luttinger and Kohn (1955) the envelope functions are limited to wave vectors inside the first Brillouin zone. Hence, the use of a local differential effective-mass equation, which gives rise to nonvanishing $`F(k)`$ outside the Brillouin zone, generates an error in the kinetic energy of order

$$(\overline{k}a)^3\left(\frac{\mathrm{}^2\overline{k}^2}{2m^{}}\right)(\overline{k}a)^3\mathrm{\Delta }\overline{V}.$$
(8)
Therefore, in the local approximation, the accuracy is limited by Eq. (8), and contributions beyond the level of (4) and (6) should for consistency be omitted.

Of course, this is not a fundamental limitation of the envelope-function method, as one can always choose to work in $`𝐤`$ space (which is a common choice for numerical work Baraff and Gershoni (1991); Gershoni et al. (1993); Winkler and Rössler (1993)). However, there is another similar source of error that arises from anticrossings of the bands explicitly included in the envelope-function model with those treated as remote perturbations. This occurs, for example, in the $`\mathrm{\Gamma }_6`$ conduction band of GaAs at about a third of the distance to the Brillouin zone boundary in the $`\mathrm{\Delta }`$ direction. Beyond such an anticrossing the model is no longer valid and an error of order (8) arises even in a $`𝐤`$-space formalism. To eliminate such errors one must enlarge the Hamiltonian by treating these remote bands explicitly, thus obtaining a full-zone $`𝐤𝐩`$ model. Cardona and Pollak (1966); Pollak et al. (1966); Wang and Zunger (1996); Mailhiot and Smith (1990)

Such a multi-band envelope-function model is also needed for medium-gap ($`\mathrm{\Delta }\overline{V}E_\text{g}^{(0)}`$) and narrow-gap ($`E_\text{g}^{(0)}\mathrm{\Delta }\overline{V}`$) systems, where $`E_\text{g}^{(0)}`$ is the energy gap of the reference crystal. The perturbative approach based on Eq. (1) can still be used, provided that $`\overline{E}_\text{g}`$ refers to the energy gap between the bands of interest and those treated as remote perturbations. Takhtamirov and Volkov Takhtamirov and Volkov (2000b) have considered the extreme narrow-gap limit for the case in which the dispersion is dominated by linear-$`k`$ terms for all energy ranges of interest, and showed that in this case the error generated by the local approximation is of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$. In such cases the terms (4) and (6) should be omitted for consistency. However, in practical situations one is more likely to encounter cases in which $`\mathrm{\Delta }\overline{V}E_\text{g}^{(0)}`$ and the $`k^2`$ dispersion terms are comparable to the linear-$`k`$ terms. Thus, in this paper all terms of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$ are retained, while those at the level of the local approximation (8) are omitted.

It should be noted that a fully self-consistent perturbation scheme is not always desirable. It would of course be impossible to achieve an accurate numerical prediction for the fine structure of the hydrogen atom without including the relativistic-mass correction and the Darwin term in addition to the spin-orbit splitting. However, for qualitative considerations one is often interested primarily in symmetry-breaking effects, in which case the former two contributions may justifiably be omitted. Likewise, although the present work retains all terms of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$, certain of these (such as the $`k^4`$ bulk dispersion terms) may possibly be omitted for applications in which the primary focus is on symmetry-breaking interface effects. \[As an example, Takhtamirov and Volkov have proposed a modelTakhtamirov and Volkov (2000a) based on certain terms of order $`(\overline{k}a)^3\mathrm{\Delta }\overline{V}`$, while neglecting larger terms of order $`(\overline{k}a)\mathrm{\Delta }\overline{V}`$ and $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$.\] However, the neglect of these terms can only be justified for specific individual applications and should not be presumed to hold in general.

The discussion here has focused on wide quantum wells, but this should not be taken to imply that effective-mass theory is inapplicable in other cases. For example, in a narrow quantum well, Burt (1992); Takhtamirov and Volkov (1999b) an effective-mass approximation can be developed along the same lines as for shallow impurities. Luttinger and Kohn (1955); Sham (1966) In this case, however, the above estimates of the interface terms and local approximation are no longer valid; see Ref. Takhtamirov and Volkov, 1999b for further details.

The remainder of this paper uses atomic units with $`\mathrm{}=m=e=1`$.

## III Quadratic response theory

### III.1 Basic definitions

The problem of interest is the Dyson eigenvalue equation (see Ref. Foreman, 2005b and references therein)

$$\frac{1}{2}^2\psi (𝐱,\omega )+V(𝐱,𝐱^{},\omega )\psi (𝐱^{},\omega )d^3x^{}=E(\omega )\psi (𝐱,\omega ),$$
(9)
in which $`\omega `$ is a complex energy parameter, $`E`$ is the complex eigenvalue, and the self-consistent potential $`V`$ is the sum of a fixed norm-conserving ionic pseudopotential $`V_{\text{ion}}`$ and the nonhermitian self-energy operator $`\mathrm{\Sigma }`$:

$$V(𝐱,𝐱^{},\omega )=V_{\text{ion}}(𝐱,𝐱^{})+\mathrm{\Sigma }(𝐱,𝐱^{},\omega ).$$
(10)
For notational simplicity, the $`\omega `$ dependence will be suppressed in most of the equations that follow. In Eq. (9), $`\psi `$ and $`V`$ are spinors, with $`V`$ having the form

$$V=𝟙𝕍_{\text{sc}}+𝝈𝐕,$$
(11)
in which $`𝟙`$ is the $`2\times 2`$ spinor unit matrix, $`V_{\text{sc}}`$ is the scalar relativistic part of $`V`$, $`𝝈`$ is the Pauli matrix, and $`𝐕`$ accounts for spin-orbit coupling. Hybertsen and Louie (1986b); Surh et al. (1991); Hemstreet et al. (1993); Theurich and Hill (2001) Note that Eq. (9) incorporates all relativistic corrections of order $`Z^2\alpha ^2`$ (where $`Z`$ is the atomic number and $`\alpha `$ is the fine-structure constant), but neglects terms of order $`\alpha ^2`$, such as spin-orbit coupling outside the atomic cores. Kleinman (1980); Bachelet and Schlüter (1982)

In a heterostructure, it is convenient to partition the ionic pseudopotential as

$$V_{\text{ion}}(𝐱,𝐱^{})=V_{\text{ion}}^{(0)}(𝐱,𝐱^{})+\mathrm{\Delta }V_{\text{ion}}(𝐱,𝐱^{}),$$
(12)
where $`V_{\text{ion}}^{(0)}`$ is the ionic pseudopotential of some periodic reference crystal (which may be a virtual crystal). This is defined as

$$V_{\text{ion}}^{(0)}(𝐱,𝐱^{})=\underset{a,j,𝐑}{}f^{aj}v_{\text{ion}}^a(𝐱𝐑_j,𝐱^{}𝐑_j),$$
(13)
in which $`v_{\text{ion}}^a(𝐱,𝐱^{})`$ is the ionic pseudopotential for atomic species $`a`$, $`f^{aj}`$ is the fractional weight associated with atom not (c) $`a`$ on site $`j`$ at position $`𝝉_j`$ in the unit cell of the reference crystal, and $`𝐑_j=𝐑+𝝉_j`$, where $`𝐑`$ is a Bravais lattice vector of the reference crystal. $`f^{aj}`$ must satisfy

$$0f^{aj}1,\underset{a}{}f^{aj}=1,$$
(14)
although the former constraint need not be strictly enforced.

The term $`\mathrm{\Delta }V_{\text{ion}}`$ is the perturbation due to the heterostructure:

$$\mathrm{\Delta }V_{\text{ion}}(𝐱,𝐱^{})=\underset{a,j,𝐑}{}\theta _𝐑^{aj}v_{\text{ion}}^a(𝐱𝐑_j,𝐱^{}𝐑_j).$$
(15)
Here $`\theta _𝐑^{aj}`$ is the change in fractional weight (relative to the reference crystal) of atom $`a`$ at position $`𝐑_j`$ in the heterostructure, which must satisfy

$$0f^{aj}+\theta _𝐑^{aj}1,\underset{a}{}\theta _𝐑^{aj}=0.$$
(16)
The constraint (16) permits one to rewrite $`\mathrm{\Delta }V_{\text{ion}}`$ as

$$\mathrm{\Delta }V_{\text{ion}}(𝐱,𝐱^{})=\underset{\alpha ,𝐑}{^{}}\theta _𝐑^\alpha \mathrm{\Delta }v_{\text{ion}}^\alpha (𝐱𝐑_\alpha ,𝐱^{}𝐑_\alpha ),$$
(17)
in which

$$\mathrm{\Delta }v_{\text{ion}}^\alpha (𝐱,𝐱^{})v_{\text{ion}}^a(𝐱,𝐱^{})v_{\text{ion}}^{\overline{a}_j}(𝐱,𝐱^{}).$$
(18)
Here $`\alpha =(a,j)`$ is a composite index, $`\overline{a}_j`$ is the label of some given atom on site $`j`$, and the prime on the summation symbol indicates that the values $`\overline{\alpha }=(\overline{a}_j,j)`$ are excluded. To simplify the interpretation of (17), the atom $`\overline{a}_j`$ is chosen to be the same as the (virtual) atom on site $`j`$ in the reference crystal.

In this paper the self-consistent potential $`V(𝐱,𝐱^{})`$ is treated using nonlinear response theory. Foreman (2005b) The fundamental assumption is that $`V`$ can be expressed as a power series in the variables $`\theta _𝐑^\alpha `$:

$$V(𝐱,𝐱^{})=V^{(0)}(𝐱,𝐱^{})+V^{(1)}(𝐱,𝐱^{})+V^{(2)}(𝐱,𝐱^{})+\mathrm{}.$$
(19)
Here $`V^{(0)}(𝐱,𝐱^{})`$ is the self-consistent potential of the reference crystal, which has the same periodicity as the ionic potential (13):

$$V^{(0)}(𝐱,𝐱^{})=V^{(0)}(𝐱+𝐑,𝐱^{}+𝐑).$$
(20)
The linear response to the heterostructure perturbation has the form

$$V^{(1)}(𝐱,𝐱^{})=\underset{\alpha ,𝐑}{^{}}\theta _𝐑^\alpha \mathrm{\Delta }v_𝐑^\alpha (𝐱,𝐱^{}),$$
(21)
while the quadratic response is

$$V^{(2)}(𝐱,𝐱^{})=\underset{\alpha ,𝐑}{^{}}\underset{\alpha ^{},𝐑^{}}{^{}}\theta _𝐑^\alpha \theta _𝐑^{}^\alpha ^{}\mathrm{\Delta }v_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐱,𝐱^{}).$$
(22)
Here the expansion coefficients in the power series are defined by

$$\begin{array}{cc}\hfill \mathrm{\Delta }v_𝐑^\alpha (𝐱,𝐱^{})& =\frac{V(𝐱,𝐱^{})}{\theta _𝐑^\alpha },\hfill \\ \hfill \mathrm{\Delta }v_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐱,𝐱^{})& =\frac{1}{2}\frac{^2V(𝐱,𝐱^{})}{\theta _𝐑^\alpha \theta _𝐑^{}^\alpha ^{}},\hfill \end{array}$$
(23)
where the derivatives are evaluated with respect to the reference crystal. These derivatives can be evaluated numerically by applying perturbations (17) with small values of $`\theta _𝐑^\alpha `$.

The leading source of error in quadratic response theory lies in the neglected cubic response term $`V^{(3)}`$. This error can be reduced somewhat by a suitable choice of reference potential. For example, in a GaAs/AlAs heterostructure, choosing $`\text{Al}_{0.5}\text{Ga}_{0.5}\text{As}`$ as a reference crystal instead of GaAs would reduce the cubic error in AlAs by a factor of 8. However, it would also create a comparable error in GaAs. Thus, the cubic error in the interface Hamiltonian for an $`\text{Al}_{0.5}\text{Ga}_{0.5}\text{As}`$ reference potential would be about $`\frac{1}{4}`$ that of a GaAs or AlAs reference potential. This is certainly an improvement, but since the quadratic response is already quite small (see Sec. III.2), it is unlikely to make much practical difference.

In the momentum representation, the Dyson equation (9) has the form

$$\frac{1}{2}k^2\psi (𝐤)+\underset{𝐤^{}}{}V(𝐤,𝐤^{})\psi (𝐤^{})=E\psi (𝐤).$$
(24)
Making use of the translation symmetry of the reference crystal, one can write the Fourier transforms of the potentials (23) as

$$\begin{array}{cc}\hfill \mathrm{\Delta }v_𝐑^\alpha (𝐤,𝐤^{})& e^{i(𝐤𝐤^{})𝐑_\alpha }\mathrm{\Delta }v^\alpha (𝐤,𝐤^{}),\hfill \\ \hfill \mathrm{\Delta }v_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐤,𝐤^{})& e^{i(𝐤𝐤^{})𝐑_\alpha }\mathrm{\Delta }v^{\alpha \alpha ^{}𝐑^{\prime \prime }}(𝐤,𝐤^{}),\hfill \end{array}$$
(25)
in which $`𝐑^{\prime \prime }=𝐑^{}𝐑`$, and the coordinate origin is $`𝐑_\alpha `$ for the modified functions on the right-hand side. The Fourier transform of the linear response (21) can therefore be written as

$$V^{(1)}(𝐤,𝐤^{})=N\underset{\alpha }{^{}}\theta ^\alpha (𝐤𝐤^{})\mathrm{\Delta }v^\alpha (𝐤,𝐤^{}).$$
(26)
Here $`N=\mathrm{\Omega }/\mathrm{\Omega }_0`$ is the number of unit cells (of volume $`\mathrm{\Omega }_0`$) in the reference crystal (of volume $`\mathrm{\Omega }`$), $`\theta ^\alpha (𝐤)`$ is the Fourier transform of $`\theta _𝐑^\alpha `$:

$`\theta _𝐑^\alpha `$ $`={\displaystyle \underset{𝐤\mathrm{\Omega }_0^{}}{}}\theta ^\alpha (𝐤)e^{i𝐤𝐑_\alpha },`$ (27a)
$`\theta ^\alpha (𝐤)`$ $`={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐑\mathrm{\Omega }}{}}\theta _𝐑^\alpha e^{i𝐤𝐑_\alpha },`$ (27b)
and $`\mathrm{\Omega }_0^{}=(2\pi )^3/\mathrm{\Omega }_0`$ is the volume of a unit cell in the reciprocal lattice. Note from (27b) that $`\theta ^\alpha (𝐤)`$ is quasi-periodic:

$$\theta ^\alpha (𝐤+𝐆)=\theta ^\alpha (𝐤)e^{i𝐆𝝉_\alpha },$$
(28)
where $`𝐆`$ is a reciprocal lattice vector of the reference crystal.

In a similar fashion, inserting (25) into the Fourier transform of (22) and relabeling $`𝐑^{}`$ gives

$$V^{(2)}(𝐤,𝐤^{})=N\underset{\alpha ,\alpha ^{},𝐑^{}}{^{}}\theta ^{\alpha \alpha ^{}𝐑^{}}(𝐤𝐤^{})\mathrm{\Delta }v^{\alpha \alpha ^{}𝐑^{}}(𝐤,𝐤^{}),$$
(29)
in which $`\theta ^{\alpha \alpha ^{}𝐑^{}}(𝐤)`$ is the Fourier transform of the pair distribution function $`\theta ^{\alpha \alpha ^{}𝐑^{}}(𝐑)\theta _𝐑^\alpha \theta _{𝐑+𝐑^{}}^\alpha ^{}`$.

### III.2 Estimation of magnitude

A crude estimate of the relative magnitudes of the linear and quadratic response can be obtained from a simple nonlinear Thomas-Fermi model, Hohenberg and Kohn (1964) with the result For

$$V^{(1)}\mathrm{\Delta }\overline{V},\frac{V^{(2)}}{V^{(1)}}\frac{\mathrm{\Delta }\overline{V}}{4ϵ_\text{F}},$$
(30)
where $`\mathrm{\Delta }\overline{V}`$ is the typical band offset defined in Sec. II, and $`ϵ_\text{F}`$ is the Fermi energy. Now the difference in screened pseudopotentials between typical III-V semiconductors is roughly $`\mathrm{\Delta }\overline{V}0.02`$–0.05 Ry, Cohen and Heine (1970) which yields the estimate $`|V^{(2)}/V^{(1)}|0.01`$. This suggests that the quadratic response is indeed very small.

This estimate is supported by the LDA calculations of Wang and Zunger Wang and Zunger (1997) for GaAs/AlAs heterostructures, in which they found (see Table I Wang and Zunger (1997)) that the interface band-mixing terms arising from the linear response were on average about 1000 times larger than those arising from the quadratic response (with the ratio between the smallest linear and largest quadratic terms being about 100, in agreement with the estimate obtained above). However, the LDA calculations for GaAs/AlAs and In0.53Ga0.47As/InP presented in Refs. Baroni et al., 1989 and Peressi et al., 1990 indicate that the quadratic density response is only about 10 times smaller than the linear density response (the linear and quadratic potentials were not given in these papers).

Nevertheless, this is still a sufficiently large ratio to establish the validity of the quadratic approximation used here. For the purposes of the perturbation scheme of Sec. II, the factor $`4ϵ_\text{F}`$ in the denominator of Eq. (30) will be treated formally as of order $`\overline{E}_\text{g}`$, so that $`\mathrm{\Delta }V^{(2)}`$ is considered to be of the same order \[$`(\mathrm{\Delta }\overline{V})^2/\overline{E}_\text{g}(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$\] as the smallest terms retained in Sec. II.

### III.3 Functional form

As shown in Ref. Foreman, 2005b, the Coulomb interaction gives rise to singularities in the linear and quadratic potentials $`\mathrm{\Delta }v(𝐤,𝐤^{})`$ when $`𝐤𝐤^{}`$ is equal to a reciprocal lattice vector. The explicit form of the linear potential is Foreman (2005b)

$$\begin{array}{c}\mathrm{\Delta }v_𝐑^\alpha (𝐤+𝐆,𝐤^{}+𝐆^{})=w_𝐑^\alpha (𝐤,𝐤^{};𝐆,𝐆^{})\hfill \\ \hfill +\mathrm{\Lambda }_{\mathrm{𝐆𝐆}^{}}(𝐤,𝐤^{})\phi _𝐑^\alpha (𝐤𝐤^{}),\end{array}$$
(31)
where the potential $`w_𝐑^\alpha (𝐤,𝐤^{};𝐆,𝐆^{})`$ is an analytic function of $`𝐤`$ and $`𝐤^{}`$ for wave vectors inside the first Brillouin zone of the reference crystal. For small values of $`𝐤`$ and $`𝐤^{}`$ it can be represented as a Taylor series, which is the basis for the effective-mass theory developed in the following section. The effective vertex function $`\mathrm{\Lambda }_{\mathrm{𝐆𝐆}^{}}(𝐤,𝐤^{})`$ is also an analytic function of $`𝐤`$ and $`𝐤^{}`$.

The singular contributions come from the screened potential $`\phi `$, which is a spin scalar of the form Foreman (2005b); not (d)

$$\phi _𝐑^\alpha (𝐪)=\frac{v_c(𝐪)n_𝐑^\alpha (𝐪)}{ϵ(𝐪)},$$
(32)
where $`ϵ(𝐪)`$ is the static electronic dielectric function of the reference crystal and

$$v_c(𝐪)=\{\begin{array}{cc}4\pi /q^2\hfill & \text{if }q0\text{,}\hfill \\ 0\hfill & \text{if }q=0\text{.}\hfill \end{array}$$
(33)
The function $`n_𝐑^\alpha (𝐪)`$ is an effective electron density containing partial contributions from the bare ionic pseudocharge and the screening charge. This is an analytic function of $`𝐪`$ (for small $`q`$) with the symmetry of site $`𝐑_\alpha `$ in the reference crystal. For isovalent substitutions in zinc-blende crystals, the leading nonanalytic terms in Eq. (32) are octopole and hexadecapole potentials proportional to $`q_xq_yq_z/q^2`$ and $`(q_x^4+q_y^4+q_z^4)/q^2`$. Foreman (2005b)

The quadratic potential has a similar form:

$$\begin{array}{c}\mathrm{\Delta }v_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐤+𝐆,𝐤^{}+𝐆^{})=w_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐤,𝐤^{};𝐆,𝐆^{})\hfill \\ \hfill +\mathrm{\Lambda }_{\mathrm{𝐆𝐆}^{}}(𝐤,𝐤^{})\phi _{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐤𝐤^{}),\end{array}$$
(34)
in which $`\mathrm{\Lambda }`$ is the same vertex function as above. In this case, however, the potential $`w_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐤,𝐤^{};𝐆,𝐆^{})`$ is not an analytic function of $`𝐤`$ and $`𝐤^{}`$. Nevertheless, to within the accuracy required here (namely, zeroth order in $`𝐤`$ and $`𝐤^{}`$), the nonanalytic part can be neglected for neutral perturbations. Note that the analytic part of the quadratic potential is analytic over a smaller region in $`𝐤`$ space than the linear potential. Whereas $`w_𝐑^\alpha (𝐤,𝐤^{};𝐆,𝐆^{})`$ was analytic for $`𝐤`$ and $`𝐤^{}`$ inside the first Brillouin zone, the analytic part of $`w_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐤,𝐤^{};𝐆,𝐆^{})`$ is analytic only over the inner “half” of the Brillouin zone (i.e., over the Brillouin zone of a crystal whose lattice constants are double those of the reference crystal).

The quadratic screened potential is defined by

$$\phi _{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐪)=\frac{v_c(𝐪)n_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐪)}{ϵ(𝐪)},$$
(35)
in which the effective density $`n_{\mathrm{𝐑𝐑}^{}}^{\alpha \alpha ^{}}(𝐪)`$ is likewise not analytic, but can be approximated as such. The leading terms here are dipole and quadrupole potentials (since the monopole term vanishes for an insulator at zero temperature). In zinc-blende crystals, only the interface dipole term is non-negligible under the approximation scheme used in this paper.

For use in Eqs. (26) and (29), one requires also expressions for the modified linear and quadratic potentials defined in Eq. (25). These are given by

$$\begin{array}{c}\mathrm{\Delta }v^\alpha (𝐤+𝐆,𝐤^{}+𝐆^{})=w_{\mathrm{𝐆𝐆}^{}}^\alpha (𝐤,𝐤^{})\hfill \\ \hfill +e^{i(𝐆𝐆^{})𝝉_\alpha }\mathrm{\Lambda }_{\mathrm{𝐆𝐆}^{}}(𝐤,𝐤^{})\phi ^\alpha (𝐤𝐤^{}),\end{array}$$
(36a)
$$\begin{array}{c}\mathrm{\Delta }v^{\alpha \alpha ^{}𝐑^{}}(𝐤+𝐆,𝐤^{}+𝐆^{})=w_{\mathrm{𝐆𝐆}^{}}^{\alpha \alpha ^{}𝐑^{}}(𝐤,𝐤^{})\hfill \\ \hfill +e^{i(𝐆𝐆^{})𝝉_\alpha }\mathrm{\Lambda }_{\mathrm{𝐆𝐆}^{}}(𝐤,𝐤^{})\phi ^{\alpha \alpha ^{}𝐑^{}}(𝐤𝐤^{}),\end{array}$$
(36b)
in which the terms on the right-hand side are defined by expressions similar to (25).

## IV Envelope-function equations

In this section, the Dyson equation (24) for a heterostructure is written in the Luttinger-Kohn representation, Luttinger and Kohn (1955) in which the basis functions are defined to be plane waves multiplied by the zone-center Bloch functions $`U_n`$ of the periodic reference crystal. Leibler (1975); Burt (1988) Here the set $`\{U_n\}`$ is chosen to be a complete set of solutions to the equation \[cf. Eq. (9)\]

$$\begin{array}{c}\frac{1}{2}^2U_n(𝐱,\omega )+V^{(\text{r},0)}(𝐱,𝐱^{},\omega )U_n(𝐱^{},\omega )d^3x^{}\hfill \\ \hfill =E_n(\omega )U_n(𝐱,\omega )\end{array}$$
(37)
that satisfy the periodic boundary conditions

$$U_n(𝐱,\omega )=U_n(𝐱+𝐑,\omega ),$$
(38)
in which $`V^{(\text{r},0)}`$ is the hermitian part of the reference potential $`V^{(0)}`$. For a general operator $`A`$, the “real” (hermitian) and “imaginary” (antihermitian) parts are defined here as

$$A^{\text{(r)}}=\frac{1}{2}(A^{}+A),A^{\text{(i)}}=\frac{i}{2}(A^{}A).$$
(39)
The solutions to Eq. (37) form a complete orthonormal set of periodic functions Luttinger and Kohn (1955) for any value of $`\omega `$. In the Fourier series representation

$$U_n(𝐱,\omega )=\underset{𝐆}{}U_{n𝐆}(\omega )e^{i𝐆𝐱},$$
(40)
the orthogonality and completeness relations are

$`{\displaystyle \underset{𝐆}{}}U_{n𝐆}^{}(\omega )U_{n^{}𝐆}(\omega )`$ $`=\delta _{nn^{}},`$ (41a)
$`{\displaystyle \underset{n}{}}U_{n𝐆}(\omega )U_{n𝐆^{}}^{}(\omega )`$ $`=𝟙\delta _{\mathrm{𝐆𝐆}^{}},`$ (41b)
in which $`U_{n𝐆}^{}`$ denotes the hermitian conjugate of the spinor $`U_{n𝐆}`$.

Note that the Bravais lattice chosen for the periodic boundary conditions in Eq. (38) need not be the same as that in (20); for certain applications it may be preferable to impose periodicity with respect to some (mathematically defined) supercell instead. For example, in treating intervalley $`\mathrm{\Gamma }`$$`X`$ coupling in semiconductors with the zinc-blende structure, it is convenient Burt (1992); Foreman (1998) to choose a nonprimitive simple cubic unit cell of volume $`\mathrm{\Omega }_0=a^3`$ (where $`a`$ is the conventional cubic lattice constant), which encompasses four primitive fcc cells. not (e) This folds the $`X`$ valleys onto the Brillouin zone center of the supercell, thereby permitting intervalley $`\mathrm{\Gamma }`$$`X`$ mixing to be described in the same notation as that for ordinary $`\mathrm{\Gamma }`$ states (although the $`\mathrm{\Gamma }`$ and $`X`$ states are of course not coupled by the $`𝐤𝐩`$ interaction).

The wave function in the Luttinger-Kohn representation (or envelope function) $`F_n(𝐤)`$ is defined as

$`F_n(𝐤)={\displaystyle \underset{𝐆}{}}U_{n𝐆}^{}\psi (𝐤+𝐆),`$ (42a)
$`\psi (𝐤+𝐆)={\displaystyle \underset{n}{}}F_n(𝐤)U_{n𝐆}.`$ (42b)
This definition is valid for any value of $`𝐤`$, but since only those values from one unit cell $`\mathrm{\Omega }_0^{}`$ are needed to determine $`\psi (𝐱)`$, it is convenient to set $`F_n(𝐤)0`$ when $`𝐤\mathrm{\Omega }_0^{}`$. Burt (1988); not (f) The Fourier transform of (42b) is then the usual exact envelope-function expansion Luttinger and Kohn (1955); Burt (1988)

$$\psi (𝐱)=\underset{n}{}F_n(𝐱)U_n(𝐱).$$
(43)
The Dyson equation in the Luttinger-Kohn representation is given by Eqs. (24), (41), and (42) as

$$\begin{array}{c}E_nF_n(𝐤)+\underset{n^{}}{}L_{nn^{}}(𝐤)F_n^{}(𝐤)\hfill \\ \hfill +\underset{n^{}}{}\underset{𝐤^{}\mathrm{\Omega }_0^{}}{}\mathrm{\Delta }V_{nn^{}}(𝐤,𝐤^{})F_n^{}(𝐤^{})=EF_n(𝐤).\end{array}$$
(44)
Here $`\mathrm{\Delta }V=VV^{(0)}`$ is the perturbation due to the heterostructure, the matrix elements of which are

$$\mathrm{\Delta }V_{nn^{}}(𝐤,𝐤^{})=\underset{𝐆,𝐆^{}}{}U_{n𝐆}^{}\mathrm{\Delta }V(𝐤+𝐆,𝐤^{}+𝐆^{})U_{n^{}𝐆^{}}.$$
(45)
The term $`L_{nn^{}}(𝐤)`$ groups together all contributions from the bulk reference crystal Hamiltonian except $`E_n`$; i.e.,

$$L_{nn^{}}(𝐤)=V_{nn^{}}^{(0)}(𝐤)V_{nn^{}}^{(\text{r},0)}(\mathrm{𝟎})+𝐤𝐩_{nn^{}}+\frac{1}{2}k^2\delta _{nn^{}},$$
(46)
in which $`V_{nn^{}}^{(0)}(𝐤)=V_{nn^{}}^{(0)}(𝐤,𝐤)`$ is the potential energy of the reference crystal and

$$𝐩_{nn^{}}=\underset{𝐆}{}𝐆U_{n𝐆}^{}U_{n^{}𝐆}$$
(47)
is the momentum matrix of the reference crystal.

Within quadratic response theory, the perturbation $`\mathrm{\Delta }V_{nn^{}}(𝐤,𝐤^{})`$ is obtained by substituting Eqs. (26), (29), and (36) into Eq. (45). The result is

$$\mathrm{\Delta }V_{nn^{}}(𝐤,𝐤^{})=W_{nn^{}}(𝐤,𝐤^{})+\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})\phi (𝐤𝐤^{}),$$
(48)
where the vertex function

$$\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})=\underset{𝐆,𝐆^{}}{}U_{n𝐆}^{}\mathrm{\Lambda }_{\mathrm{𝐆𝐆}^{}}(𝐤,𝐤^{})U_{n^{}𝐆^{}}$$
(49)
is again an analytic function of $`𝐤`$ and $`𝐤^{}`$. The $`W`$ term is defined by $`W=W^{(1)}+W^{(2)}`$, where

$$W_{nn^{}}^{(1)}(𝐤,𝐤^{})=\underset{\alpha }{^{}}\theta ^\alpha (𝐤𝐤^{})W_{nn^{}}^\alpha (𝐤,𝐤^{}),$$
(50)
$$W_{nn^{}}^\alpha (𝐤,𝐤^{})=N\underset{𝐆,𝐆^{}}{}U_{n𝐆}^{}w_{\mathrm{𝐆𝐆}^{}}^\alpha (𝐤,𝐤^{})U_{n^{}𝐆^{}}e^{i(𝐆^{}𝐆)𝝉_\alpha }.$$
(51)
The screened potential $`\phi `$ is defined by

$$\phi (𝐪)=\frac{v_c(𝐪)n(𝐪)}{ϵ(𝐪)},$$
(52)
where $`n=n^{(1)}+n^{(2)}`$ is an effective electron density for the heterostructure perturbation:

$$n^{(1)}(𝐪)=N\underset{\alpha }{^{}}\theta ^\alpha (𝐪)n^\alpha (𝐪).$$
(53)
The quadratic contributions $`W^{(2)}`$ and $`n^{(2)}`$ are given by obvious generalizations \[see Eq. (36)\] of the above expressions. The physical interpretation of these results is considered below.

## V Power series expansions

In this section, power series expansions are used to obtain approximate expressions for the Hamiltonian matrix elements in the envelope-function equations (44). This approximation is justified by the existence Burt (1992, 1994); Xia (1989); Wang et al. (1997); Wang and Zunger (1997, 1999) of slowly varying envelope functions $`F_n(𝐱)`$, for which $`F_n(𝐤)`$ is negligible unless $`𝐤`$ is small. This expansion provides a starting point for the development of an approximate effective-mass theory, and also assists in the physical interpretation of the various terms in the Hamiltonian.

### V.1 Reference crystal Hamiltonian

The leading terms in the bulk Hamiltonian $`L_{nn^{}}(𝐤)`$ of Eq. (46) are derived from a Taylor series expansion of $`V_{nn^{}}^{(0)}(𝐤)`$, and are given through terms of the fourth order in $`k`$ by

$$\begin{array}{c}L_{nn^{}}(𝐤)=k_\lambda \pi _{nn^{}}^\lambda +k_\lambda k_\mu \stackrel{~}{D}_{nn^{}}^{\lambda \mu }+k_\lambda k_\mu k_\kappa \stackrel{~}{C}_{nn^{}}^{\lambda \mu \kappa }\hfill \\ \hfill +k_\lambda k_\mu k_\kappa k_\nu \stackrel{~}{Q}_{nn^{}}^{\lambda \mu \kappa \nu }.\end{array}$$
(54)
Here a summation with respect to the Cartesian indices $`\lambda ,\mu ,\kappa ,\text{and }\nu `$ is implicit, and the coefficients are

$$\begin{array}{cc}\hfill \pi _{nn^{}}^\lambda & =p_{nn^{}}^\lambda +\left(\frac{V_{nn^{}}^{(\text{r},0)}(𝐤)}{k_\lambda }\right)_{𝐤=\mathrm{𝟎}},\hfill \\ \hfill \stackrel{~}{D}_{nn^{}}^{\lambda \mu }& =\frac{1}{2}\left(\delta _{\lambda \mu }\delta _{nn^{}}+\frac{^2V_{nn^{}}^{(\text{r},0)}(𝐤)}{k_\lambda k_\mu }\right)_{𝐤=\mathrm{𝟎}},\hfill \\ \hfill \stackrel{~}{C}_{nn^{}}^{\lambda \mu \kappa }& =\frac{1}{3!}\left(\frac{^3V_{nn^{}}^{(\text{r},0)}(𝐤)}{k_\lambda k_\mu k_\kappa }\right)_{𝐤=\mathrm{𝟎}},\hfill \\ \hfill \stackrel{~}{Q}_{nn^{}}^{\lambda \mu \kappa \nu }& =\frac{1}{4!}\left(\frac{^4V_{nn^{}}^{(\text{r},0)}(𝐤)}{k_\lambda k_\mu k_\kappa k_\nu }\right)_{𝐤=\mathrm{𝟎}}.\hfill \end{array}$$
(55)
The various derivatives of $`V_{nn^{}}^{(\text{r},0)}(𝐤)`$ account for contributions from the nonlocal part of the potential energy to the dispersion relation of the reference crystal. The quantity $`\pi _{nn^{}}^\lambda `$ is the kinetic momentum matrix of the reference crystal, whereas the other terms give partial contributions (see Sec. VI for the remaining contributions) to the effective-mass tensor and the cubic and quartic dispersion terms of the reference crystal.

In principle, $`L_{nn^{}}(𝐤)`$ should include contributions from the antihermitian part of the self-energy $`\mathrm{\Sigma }_{nn^{}}^{\text{(i,0)}}(𝐤,\omega )=V_{nn^{}}^{\text{(i,0)}}(𝐤,\omega )`$. However, in Appendix A it is shown that, for energies $`\omega `$ near the band gap of the reference crystal, the contributions from $`\mathrm{\Sigma }^{\text{(i)}}`$ are much smaller than the smallest terms retained in the present approximation scheme. Such contributions were therefore neglected in the above expressions, and are likewise neglected in subsequent analysis of the heterostructure perturbation. However, in any calculation where it is desired to include the effects of a finite quasiparticle lifetime, the dominant terms may be restored to leading order by replacing $`E_n`$ with

$$E_n(\omega )E_n(\omega )+i\mathrm{\Sigma }_{nn}^{(\text{i},0)}(𝐤=\mathrm{𝟎},\omega ),$$
(56)
and then retaining the imaginary part only to first order in perturbation theory in all subsequent analysis.

### V.2 Linear heterostructure potential

A similar Taylor series expansion technique is useful for the terms $`W_{nn^{}}(𝐤,𝐤^{})`$ and $`\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})`$ in the heterostructure perturbation (48). This subsection begins by considering the simple special case in which the screened atomic pseudopotentials $`w_{\mathrm{𝐆𝐆}^{}}^\alpha (𝐤,𝐤^{})`$ in Eqs. (36a) and (51) have the form of a local potential; i.e.,

$$w_{\mathrm{𝐆𝐆}^{}}^\alpha (𝐤,𝐤^{})=w^\alpha (𝐤𝐤^{}+𝐆𝐆^{}).$$
(57)
Such would be the case, for example, in an LDA calculation based on local ionic pseudopotentials. In this case, $`W_{nn^{}}^{(1)}(𝐤,𝐤^{})`$ is also a local potential of the form $`W_{nn^{}}^{(1)}(𝐤,𝐤^{})=W_{nn^{}}^{(1)}(𝐤𝐤^{})`$. This simplification makes it easier to grasp the physical significance of the power series expansion, and also facilitates a comparison between the present theory and earlier envelope-function models based on local empirical pseudopotentials.

#### V.2.1 Local analytic terms

The approximation technique used here for $`W_{nn^{}}^{(1)}(𝐤𝐤^{})`$ is to expand the term $`W_{nn^{}}^\alpha (𝐤𝐤^{})`$ on the right-hand side of Eq. (50) in a Taylor series in $`𝐪𝐤𝐤^{}`$:

$$W_{nn^{}}^{(1)}(𝐪)=\underset{\alpha }{^{}}\theta ^\alpha (𝐪)[W_{nn^{}}^\alpha +iq_\lambda (Z_{nn^{}}^\lambda )^\alpha q_\lambda q_\mu (Y_{nn^{}}^{\lambda \mu })^\alpha ].$$
(58)
Here the terms $`W_{nn^{}}^\alpha `$, $`(Z_{nn^{}}^\lambda )^\alpha `$, and $`(Y_{nn^{}}^{\lambda \mu })^\alpha `$ are $`𝐪`$-independent expansion coefficients. The series (58) has been truncated at the second order in $`q`$ because such terms are of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$, which are the smallest corrections permitted in the present perturbation scheme.

The result (58) may be used directly in a $`𝐤`$-space envelope-function calculation based on Eq. (44), in which the $`𝐤`$ values are expressly limited to the unit cell $`\mathrm{\Omega }_0^{}`$. However, to obtain a local differential equation, one must allow $`𝐤`$ to range over all possible values. Since $`\theta ^\alpha (𝐤)`$ is quasi-periodic \[see Eq. (28)\], this local approximation will generate large-$`𝐤`$ terms in the envelope functions unless a $`𝐤`$-space cutoff is introduced. If this is done, the Fourier transform of (58) is the local potential

$$W_{nn^{}}^{(1)}(𝐱)=\underset{\alpha }{^{}}[W_{nn^{}}^\alpha +(Z_{nn^{}}^\lambda )^\alpha _\lambda +(Y_{nn^{}}^{\lambda \mu })^\alpha _\lambda _\mu ]\theta ^\alpha (𝐱),$$
(59)
in which $`_\lambda =/x_\lambda `$ and \[cf. Eq. (27a)\]

$$\theta ^\alpha (𝐱)=\underset{𝐤}{}B(𝐤)\theta ^\alpha (𝐤)e^{i𝐤𝐱},$$
(60)
where the cutoff function $`B(𝐤)`$ is defined in Appendix B. From Eq. (16), the constraint

$$\underset{a}{}\theta ^{aj}(𝐱)=0$$
(61)
is satisfied for any choice of $`B(𝐤)`$.

The physical significance of the result (59) can be appreciated by considering a specific example such as a (001) GaAs/AlAs heterojunction. In this case, as will be shown in Sec. VII, the function $`\theta ^\alpha (𝐱)`$ depends only on the $`z`$ coordinate and behaves like a smooth step function at the interface. The spatial derivatives in Eq. (59) therefore generate finite-width $`\delta `$-like terms at the interface, with the $`Z`$ term proportional to $`\delta (z)`$ and the $`Y`$ term proportional to $`\delta ^{}(z)`$. Hence, Eq. (59) provides a first example of the interface band-mixing terms alluded to previously in Eqs. (5) and (6).

The physical origin of these terms can be understood by going back one step further in the derivation. From Eq. (51), it is clear that the Taylor series expansion of $`W_{nn^{}}^\alpha (𝐪)`$ in Eq. (58) is equivalent to a Taylor series expansion of the screened atomic pseudopotential (57) with respect to $`𝐪=𝐤𝐤^{}`$. Hence, the physical origin of the linear and quadratic (in $`q`$) terms in Eq. (58) is simply the finite slope and curvature of the screened atomic pseudopotentials in momentum space.

This demonstrates that the Hamiltonian of a heterostructure depends not just on the values of the atomic pseudopotentials $`w^\alpha (𝐤)`$ at the reciprocal lattice vectors $`𝐆`$, but also at a range of $`𝐤`$ values in a finite neighborhood of each $`𝐆`$ (the size of the neighborhood depending on how rapidly varying the envelope function is). The necessity for an accurate fitting of empirical pseudopotentials over a range of $`𝐤`$ values has been emphasized particularly in the work of Mäder and Zunger. Mäder and Zunger (1994) The truncated expansion (58) shows that the present perturbation scheme relies for its accuracy upon the validity of a quadratic extrapolation of $`w^\alpha (𝐤)`$ in the neighborhood of each $`𝐆`$. Inspection of the form of typical atomic pseudopotentials Cohen and Heine (1970); Cohen and Chelikowsky (1989) shows this to be a good approximation.

The definition (58) of the interface band-mixing parameters $`Z`$ and $`Y`$ highlights a significant difference between the present approach and previous derivations of envelope-function Hamiltonians. In previous derivations, the heterostructure potential was chosen to have the form

$$V(𝐱)=\underset{l}{}V^l(𝐱)\theta ^l(𝐱),$$
(62)
in which $`V^l(𝐱)`$ is the periodic potential for the bulk material $`l`$, and $`\theta ^l(𝐱)`$ is a form factor determining the composition of the heterostructure. Within this model, $`V(𝐱)`$ depends only upon the atomic pseudopotentials at the reciprocal lattice vectors $`𝐆`$, with the magnitude of the interface $`\delta `$-like terms determined by $`\theta ^l(𝐱)`$. In the limit of slowly varying $`\theta ^l(𝐱)`$, the interface terms vanish.

The present results show that this behavior is an unphysical artifact of the model (62). In the present theory, the strength of the interface terms is the same (in the virtual crystal approximation) for slowly graded and abrupt heterostructures, with the rate of change of $`\theta ^\alpha (𝐱)`$ affecting only the width of the interface terms.

#### V.2.2 Nonlocal analytic terms

The more general case of a nonlocal potential involves only a straightforward extension of these results. The matrix element $`W_{nn^{}}^\alpha (𝐤,𝐤^{})`$ of Eq. (51) is expanded to second order in $`𝐤`$ and $`𝐤^{}`$ as follows:

$$\begin{array}{c}W_{nn^{}}^\alpha (𝐤,𝐤^{})=W_{nn^{}}^\alpha +k_\lambda (\stackrel{~}{J}_{nn^{}}^\lambda )^\alpha +k_\lambda ^{}[(\stackrel{~}{J}_{n^{}n}^\lambda )^\alpha ]^{}\hfill \\ \hfill +k_\lambda k_\mu (\stackrel{~}{M}_{nn^{}}^{\lambda \mu })^\alpha +k_\lambda ^{}k_\mu ^{}[(\stackrel{~}{M}_{n^{}n}^{\mu \lambda })^\alpha ]^{}\\ \hfill +k_\lambda k_\mu ^{}(\stackrel{~}{R}_{nn^{}}^{\lambda \mu })^\alpha .\end{array}$$
(63)
Here the antihermitian part of the self-energy was neglected, as discussed above Eq. (56) and in Appendix A. This approximation simplifies the expansion by providing a relationship between (for example) the coefficients of $`k_\lambda `$ and $`k_\lambda ^{}`$. The physical interpretation of the various expansion coefficients is discussed below (in Sec. VI), after perturbation theory has been used to eliminate the coupling to remote bands.

When this expansion is substituted into Eq. (50) for $`W_{nn^{}}^{(1)}(𝐤,𝐤^{})`$, the result can be written as

$$\begin{array}{c}W_{nn^{}}^{(1)}(𝐤,𝐤^{})=W_{nn^{}}^{(1)}(𝐪)+k_\lambda \stackrel{~}{J}_{nn^{}}^\lambda (𝐪)+k_\lambda ^{}[\stackrel{~}{J}_{n^{}n}^\lambda (𝐪)]^{}\hfill \\ \hfill +k_\lambda k_\mu \stackrel{~}{M}_{nn^{}}^{\lambda \mu }(𝐪)+k_\lambda ^{}k_\mu ^{}[\stackrel{~}{M}_{n^{}n}^{\mu \lambda }(𝐪)]^{}\\ \hfill +k_\lambda k_\mu ^{}\stackrel{~}{R}_{nn^{}}^{\lambda \mu }(𝐪),\end{array}$$
(64)
in which $`𝐪=𝐤𝐤^{}`$ as before, and the various $`𝐪`$-dependent functions are defined by expressions of the form

$$\stackrel{~}{R}_{nn^{}}^{\lambda \mu }(𝐪)=\underset{\alpha }{^{}}\theta ^\alpha (𝐪)(\stackrel{~}{R}_{nn^{}}^{\lambda \mu })^\alpha .$$
(65)
These functions all have a step-function-like behavior in $`𝐱`$ space at a heterojunction (see Sec. VII). Since $`\theta _𝐑^\alpha `$ is real, they also have the hermiticity and symmetry properties

$$\begin{array}{cc}\hfill [W_{nn^{}}^{(1)}(𝐪)]^{}& =W_{n^{}n}^{(1)}(𝐪),\hfill \\ \hfill [\stackrel{~}{R}_{nn^{}}^{\lambda \mu }(𝐪)]^{}& =\stackrel{~}{R}_{n^{}n}^{\mu \lambda }(𝐪),\hfill \\ \hfill \stackrel{~}{M}_{nn^{}}^{\lambda \mu }(𝐪)& =\stackrel{~}{M}_{nn^{}}^{\mu \lambda }(𝐪).\hfill \end{array}$$
(66)
#### V.2.3 Nonanalytic terms

The term $`\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})\phi ^{(1)}(𝐪)`$ describing the nonanalytic contributions to Eq. (48) is handled in much the same way. The vertex function $`\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})`$ is expanded in a Taylor series of the form (63), while the linear electron density in Eq. (53) is expanded as

$$n^\alpha (𝐪)=q_\lambda n_\lambda ^\alpha +q_\lambda q_\mu n_{\lambda \mu }^\alpha +q_\lambda q_\mu q_\kappa n_{\lambda \mu \kappa }^\alpha +q_\lambda q_\mu q_\kappa q_\nu n_{\lambda \mu \kappa \nu }^\alpha ,$$
(67)
which is valid for neutral perturbations. In a heterovalent zinc-blende heterostructure described by neutral perturbations, Foreman (2005b); not (a) the dipole and traceless quadrupole terms $`n_\lambda ^\alpha `$ and $`n_{\lambda \mu }^\alpha \frac{1}{3}n_{\nu \nu }^\alpha \delta _{\lambda \mu }`$ are nonvanishing only at interfaces. For isovalent substitutions in zinc-blende, $`n_\lambda ^\alpha =0`$ and $`n_{\lambda \mu }^\alpha =\frac{1}{3}n_{\nu \nu }^\alpha \delta _{\lambda \mu }`$ everywhere, so the contribution from the latter term is analytic and can be absorbed into the definition of the analytic potential (63). The remaining terms $`n_{\lambda \mu \kappa }^\alpha `$ and $`n_{\lambda \mu \kappa \nu }^\alpha `$ describe octopole and hexadecapole moments. Foreman (2005b)

The inverse dielectric function $`ϵ^1(𝐪)`$ in Eq. (52) can also be expanded in a power series. For isovalent perturbations in zinc-blende this has no qualitative significance, since the leading correction merely renormalizes the hexadecapole term $`n_{\lambda \mu \kappa \nu }^\alpha `$. Foreman (2005b) However, the wave vector dependence of $`ϵ^1(𝐪)`$ does generate a qualitatively new contribution from the interface dipole term in Eq. (67); see Eq. (6.23) of Ref. Foreman, 2005b for the explicit form of this term.

With the expansion (67), the effective density (53) has the same form as that derived above for the local potential in Eq. (58). In coordinate space, it involves a series of derivatives of $`\theta ^\alpha (𝐱)`$, similar to the result shown in Eq. (59). This is the same as the usual multipole expansion of the macroscopic charge density in classical electromagnetism. Jackson (1999) In particular, note that for a perturbation consisting of a single impurity atom, the function $`\theta ^\alpha (𝐱)`$ is just the macroscopic average of a Dirac $`\delta `$ function, in complete agreement with the expressions given in Ref. Jackson, 1999.

Since the Taylor series expansion for $`\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})`$ is identical in form to that for $`W_{nn^{}}^\alpha (𝐤,𝐤^{})`$, the contribution from this term will not be written out explicitly here. All subsequent perturbation theory analysis for the two terms is formally identical, except that the vertex function is multiplied by an extra factor of $`\phi (𝐪)`$.

### V.3 Quadratic heterostructure potential

#### V.3.1 Analytic terms

Because the quadratic response is already of order $`(\mathrm{\Delta }\overline{V})^2/\overline{E}_\text{g}`$, the quadratic version of Eq. (50) can be replaced by the zeroth-order approximation

$`W_{nn^{}}^{(2)}(𝐤,𝐤^{})`$ $`=\underset{\alpha ,\alpha ^{},𝐑^{}}{{\displaystyle ^{}}}\theta ^{\alpha \alpha ^{}𝐑^{}}(𝐪)W_{nn^{}}^{\alpha \alpha ^{}𝐑^{}}`$
$`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪),`$ (68)
in which

$$W_{nn^{}}^{\alpha \alpha ^{}𝐑^{}}=\underset{k,k^{}0}{lim}W_{nn^{}}^{\alpha \alpha ^{}𝐑^{}}(𝐤,𝐤^{}),$$
(69)
where the limit is well defined for neutral perturbations. Foreman (2005b) This contributes a local potential-energy term similar to that given by $`W_{nn^{}}^{(1)}(𝐪)`$ in Eq. (64). At a heterojunction, the functional dependence in $`𝐱`$ space is similar to that of a smooth (macroscopically averaged) step function, with possible deviations in the vicinity of the junction.

Now since $`W_{nn^{}}^{\alpha \alpha ^{}𝐑^{}}`$ depends only on the part of $`W_{nn^{}}^{\alpha \alpha ^{}𝐑^{}}(𝐤,𝐤^{})`$ that is analytic in $`𝐤`$ and $`𝐤^{}`$, it is a short-range quantity that is significant only when $`𝐑^{}+𝝉_\alpha ^{}𝝉_\alpha `$ is comparable to the lattice constant $`a`$. Therefore, at a heterojunction, replacing $`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)`$ with an ideal step function would generate an error of order $`\overline{k}a`$ in a term of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$. Hence, the error is of the same order as the local approximation (8) and can be neglected.

Thus, within the present perturbation scheme there is no interface contribution from $`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)`$. This is an important result, as it simplifies the analysis of interface effects in later sections of this paper.

Within the context of an empirical pseudopotential model, Mäder and Zunger (1994) the nonlinear bulk term $`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)`$ derived here can be represented as an environment dependence of the screened empirical pseudopotential, in which (for example) the pseudopotential for an As atom in GaAs is different from that for an As atom in AlAs. The importance of accounting for such effects has been emphasized by Mäder and Zunger. Mäder and Zunger (1994)

#### V.3.2 Nonanalytic terms

The leading terms in the quadratic density are the dipole and quadrupole terms

$$n^{\alpha \alpha ^{}𝐑^{}}(𝐪)=q_\lambda n_\lambda ^{\alpha \alpha ^{}𝐑^{}}+q_\lambda q_\mu n_{\lambda \mu }^{\alpha \alpha ^{}𝐑^{}},$$
(70)
which again is valid for neutral perturbations. (Here the constant term vanishes even for charged perturbations as long as the system is insulating, Foreman (2005b) but in this case there is an additional nonanalytic term of order $`q^2`$.Foreman (2005b)) In zinc-blende materials, the net contributions to $`n^{(2)}(𝐪)`$ from $`n_\lambda ^{\alpha \alpha ^{}𝐑^{}}`$ and the traceless part of $`n_{\lambda \mu }^{\alpha \alpha ^{}𝐑^{}}`$ both vanish in the bulk regions of a heterostructure. The interface part of $`n_{\lambda \mu }^{\alpha \alpha ^{}𝐑^{}}`$ is negligible under the current approximation scheme, while the bulk part can be absorbed into the definition of the analytic potential $`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)`$. Therefore, only the interface dipole term remains under the current approximation scheme.

### V.4 Two-dimensional systems

In a heterostructure (such as a quantum well) with two-dimensional translation symmetry, the nonanalytic terms arising from $`\phi (𝐪)`$ have a particularly simple form. Let the dimensionless coordinates $`x_i`$ and $`k_j`$ be defined by $`𝐱=x_i𝐚_i`$ and $`𝐤=k_j𝐛_j`$, where $`𝐚_i`$ and $`𝐛_j`$ are basis vectors for the direct and reciprocal lattices of the reference crystal, with $`𝐚_i𝐛_j=2\pi \delta _{ij}`$. In these coordinates, the lattice sites are defined by $`x_i=R_i`$ and $`k_j=G_j`$, where $`R_i`$ and $`G_j`$ are integers. In such a two-dimensional system, one can choose $`𝐚_1`$ and $`𝐚_2`$ to lie parallel to the junction plane, so that the atomic distribution function $`\theta _𝐑^\alpha `$ is independent of $`R_1`$ and $`R_2`$. The Fourier transform $`\theta ^\alpha (𝐤)`$ then has the form

$$\theta ^\alpha (𝐤)=\theta ^\alpha (k_3)\underset{𝐆_{}}{}\delta _{𝐤_{}𝐆_{}}e^{i𝐆_{}𝝉_\alpha },$$
(71)
where $`𝐤_{}=k_1𝐛_1+k_2𝐛_2`$. For small $`𝐤_{}`$, the $`𝐤_{}`$ dependence is simply $`\delta _{𝐤_{}\mathrm{𝟎}}`$. The same conclusion holds for the pair distribution function $`\theta ^{\alpha \alpha ^{}𝐑^{}}(𝐤)`$.

Therefore, the nonanalytic potential (52) has the form $`\phi (𝐪)=\phi (q_3)\delta _{𝐪_{}\mathrm{𝟎}}`$, which is independent of the direction of $`𝐪`$. As a result, all terms in $`\phi (𝐪)`$ except the monopole and dipole terms reduce to analytic functions of $`q_3`$, which can be absorbed into the definition of the analytic potential $`W`$. In regard to the monopole and dipole terms, this paper considers only neutral perturbations in crystals with no bulk dipole moment. Thus, the only nonanalytic contributions are the $`1/q_3`$ terms generated by the interface dipoles in Eqs. (67) and (70).

These have the same $`q_3`$ dependence as the Fourier transform of a step function, and merely add extra terms to the band offsets at a heterojunction. Hence, the interface dipole contributions can be absorbed into the definition of $`W_{nn^{}}^{(1)}(𝐪)`$ and $`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)`$. Note that in a no-common-atom system such as InAs/GaSb, the contribution from the quadratic interface dipole to the band offset has a different value for GaAs-like and InSb-like interfaces. Dandrea et al. (1992)

In summary, for the material systems considered in this paper, the nonanalytic potential $`\phi (𝐪)`$ does not contribute anything qualitatively new in a heterostructure with two-dimensional translation symmetry. Only the interface dipole term is truly nonanalytic, and that can be absorbed into the definition of the band offsets (although this contribution does depend on the microscopic structure of the interfaceDandrea et al. (1992)). Long-range potentials arising from the direction dependence of $`\phi (𝐪)`$ appear only in structures with lower translation symmetry, such as quantum wires and dots.

## VI Elimination of interband coupling

In this section, perturbation theory Luttinger and Kohn (1955); Bir and Pikus (1974); Winkler (2003) is used to derive a multi-band effective-mass Hamiltonian from the infinite-dimensional matrix equations (44). The method is outlined briefly here; for further details see Refs. Bir and Pikus, 1974 and Winkler, 2003.

The zone-center states $`n`$ of the reference crystal are divided into a class $`A`$ containing the states of interest, and a class $`B`$ containing all other states. The total Hamiltonian is written as $`H=H_0+H^{}`$, where $`H_0`$ has matrix elements $`(H_0)_{mm^{}}=E_m\delta _{mm^{}}`$, and $`m=(n,𝐤)`$ is a composite index. The unperturbed Hamiltonian $`H_0`$ is assumed to be hermitian, but the perturbation $`H^{}`$ need not be. A similarity transformation $`\overline{H}=e^SHe^S`$ is used to eliminate the coupling between sets $`A`$ and $`B`$ to any desired order in the perturbation $`H^{}`$. (The transformation is unitary if $`H^{}`$ is hermitian, as is approximately the case here.) This yields a finite-dimensional effective Hamiltonian for states $`m,m^{}A`$, which is given explicitly (to third order in $`H^{}`$) in Eq. (113) of Appendix C.

For the case considered here, the matrix elements of $`H_0`$ and $`H^{}`$ are given by Eq. (44) as $`n𝐤|H_0|n^{}𝐤^{}=E_n\delta _{nn^{}}\delta _{\mathrm{𝐤𝐤}^{}}`$ and

$$\begin{array}{c}n𝐤|H^{}|n^{}𝐤^{}=L_{nn^{}}(𝐤)\delta _{\mathrm{𝐤𝐤}^{}}+W_{nn^{}}^{(1)}(𝐤,𝐤^{})+\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)\hfill \\ \hfill +\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})\phi (𝐪),\end{array}$$
(72)
in which $`L_{nn^{}}(𝐤)`$ is defined in Eq. (54), $`W_{nn^{}}^{(1)}(𝐤,𝐤^{})`$ in Eq. (64), and $`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)`$ in Eq. (68). Upon inserting these matrix elements into Eq. (113), one obtains the effective-mass Hamiltonian (for $`n,n^{}A`$)

$$\begin{array}{c}n𝐤|\overline{H}|n^{}𝐤^{}=(E_n\delta _{nn^{}}+k_\lambda \pi _{nn^{}}^\lambda +k_\lambda k_\mu D_{nn^{}}^{\lambda \mu }+k_\lambda k_\mu k_\kappa C_{nn^{}}^{\lambda \mu \kappa }+k_\lambda k_\mu k_\kappa k_\nu Q_{nn^{}}^{\lambda \mu \kappa \nu })\delta _{\mathrm{𝐤𝐤}^{}}+W_{nn^{}}^{(1)}(𝐪)+W_{nn^{}}^{(2)}(𝐪)\hfill \\ \hfill +k_\lambda J_{nn^{}}^\lambda (𝐪)+k_\lambda ^{}[J_{n^{}n}^\lambda (𝐪)]^{}+k_\lambda k_\mu M_{nn^{}}^{\lambda \mu }(𝐪)+k_\lambda ^{}k_\mu ^{}[M_{n^{}n}^{\mu \lambda }(𝐪)]^{}+k_\lambda k_\mu ^{}R_{nn^{}}^{\lambda \mu }(𝐪)\\ \hfill +[\mathrm{\Lambda }_{nn^{}}+k_\lambda \widehat{J}_{nn^{}}^\lambda +k_\lambda ^{}(\widehat{J}_{n^{}n}^\lambda )^{}+k_\lambda k_\mu \widehat{M}_{nn^{}}^{\lambda \mu }+k_\lambda ^{}k_\mu ^{}(\widehat{M}_{n^{}n}^{\mu \lambda })^{}+k_\lambda k_\mu ^{}\widehat{R}_{nn^{}}^{\lambda \mu }]\phi (𝐪).\end{array}$$
(73)
In this expression, $`\pi _{nn^{}}^\lambda `$ is the kinetic momentum matrix (55) of the reference crystal, $`2D_{nn^{}}^{\lambda \mu }`$ is the inverse effective mass tensor of the reference crystal \[see Eq. (114)\], and $`C_{nn^{}}^{\lambda \mu \kappa }`$ and $`Q_{nn^{}}^{\lambda \mu \kappa \nu }`$ are the coefficients of the cubic and quartic dispersion terms in the reference crystal \[see Eqs. (115) and (116)\]. These are just renormalized versions of the quantities $`\stackrel{~}{D}_{nn^{}}^{\lambda \mu }`$, $`\stackrel{~}{C}_{nn^{}}^{\lambda \mu \kappa }`$, and $`\stackrel{~}{Q}_{nn^{}}^{\lambda \mu \kappa \nu }`$ defined previously in Eq. (55).

Likewise, the functions $`W_{nn^{}}^{(2)}(𝐪)`$, $`J_{nn^{}}^\lambda (𝐪)`$, $`M_{nn^{}}^{\lambda \mu }(𝐪)`$, and $`R_{nn^{}}^{\lambda \mu }(𝐪)`$ are all renormalized versions of quantities defined previously. The renormalized functions are given explicitly in Appendix C.

The terms multiplying the screened nonanalytic potential $`\phi (𝐪)`$ are derived from the Taylor series expansion and $`𝐤𝝅`$ renormalization of the vertex function $`\mathrm{\Lambda }_{nn^{}}(𝐤,𝐤^{})`$. The various constant coefficients (e.g., $`\widehat{R}_{nn^{}}^{\lambda \mu }`$) are defined in the same way as the analogous $`𝐪`$-dependent functions \[e.g., $`R_{nn^{}}^{\lambda \mu }(𝐪)`$\] given in Appendix C, but with $`W_{nn^{}}^{(1)}(𝐪)`$ replaced by the constant $`\mathrm{\Lambda }_{nn^{}}=\mathrm{\Lambda }_{nn^{}}(\mathrm{𝟎},\mathrm{𝟎})`$.

Equation (73) can now be rearranged Leibler (1975, 1977) and Fourier transformed to obtain the effective-mass Hamiltonian

$$\begin{array}{c}\overline{H}_{nn^{}}(𝐱,𝐩)=E_n\delta _{nn^{}}+\{p_\lambda ,\pi _{nn^{}}^\lambda (𝐱)\}+\{p_\lambda p_\mu ,D_{nn^{}}^{\{\lambda \mu \}}(𝐱)\}+p_\lambda p_\mu p_\kappa C_{nn^{}}^{\lambda \mu \kappa }+p_\lambda p_\mu p_\kappa p_\nu Q_{nn^{}}^{\lambda \mu \kappa \nu }+W_{nn^{}}^{(1)}(𝐱)+W_{nn^{}}^{(2)}(𝐱)\hfill \\ \hfill +_\lambda Z_{nn^{}}^\lambda (𝐱)+_\lambda _\mu Y_{nn^{}}^{\lambda \mu }(𝐱)+\{p_\lambda ,_\mu \mathrm{\Gamma }_{nn^{}}^{\lambda \mu }(𝐱)\}+[_\mu \mathrm{\Phi }_{nn^{}}^{\lambda \mu }(𝐱)]p_\lambda \\ \hfill +\mathrm{\Lambda }_{nn^{}}\phi (𝐱)+\widehat{\pi }_{nn^{}}^\lambda \{p_\lambda ,\phi (𝐱)\}+\widehat{D}_{nn^{}}^{\{\lambda \mu \}}\{p_\lambda p_\mu ,\phi (𝐱)\}\\ \hfill +\widehat{Z}_{nn^{}}^\lambda _\lambda \phi (𝐱)+\widehat{Y}_{nn^{}}^{\lambda \mu }_\lambda _\mu \phi (𝐱)+\widehat{\mathrm{\Gamma }}_{nn^{}}^{\lambda \mu }\{p_\lambda ,_\mu \phi (𝐱)\}+\widehat{\mathrm{\Phi }}_{nn^{}}^{\lambda \mu }[_\mu \phi (𝐱)]p_\lambda ,\end{array}$$
(74)
in which $`\{A,B\}=\{AB\}=\frac{1}{2}(AB+BA)`$ is the symmetrized product, $`𝐩=i`$ is the momentum operator, and $`_\lambda =/x_\lambda `$ acts only on the function immediately to its right.

All of the $`𝐱`$-dependent functions are defined in terms of the $`𝐤`$-space cutoff (60). The quantity

$$\pi _{nn^{}}^\lambda (𝐱)=\pi _{nn^{}}^\lambda +J_{nn^{}}^\lambda (𝐱)+[J_{n^{}n}^\lambda (𝐱)]^{}$$
(75)
is the material-dependent kinetic momentum matrix for the heterostructure, while

$$D_{nn^{}}^{\lambda \mu }(𝐱)=D_{nn^{}}^{\lambda \mu }+M_{nn^{}}^{\lambda \mu }(𝐱)+[M_{n^{}n}^{\mu \lambda }(𝐱)]^{}+R_{nn^{}}^{\lambda \mu }(𝐱)$$
(76)
is half the inverse effective mass tensor for the heterostructure. This has symmetric and antisymmetric parts:

$$\begin{array}{cc}\hfill D_{nn^{}}^{\{\lambda \mu \}}(𝐱)& =\frac{1}{2}[D_{nn^{}}^{\lambda \mu }(𝐱)+D_{nn^{}}^{\mu \lambda }(𝐱)],\hfill \\ \hfill D_{nn^{}}^{[\lambda \mu ]}(𝐱)& =\frac{1}{2}[D_{nn^{}}^{\lambda \mu }(𝐱)D_{nn^{}}^{\mu \lambda }(𝐱)],\hfill \end{array}$$
(77)
although the antisymmetric part has a nonvanishing contribution only in the presence of a magnetic field.

The functions $`Z`$, $`Y`$, $`\mathrm{\Gamma }`$, and $`\mathrm{\Phi }`$ in Eq. (74) are all interface terms. The first two terms

$`Z_{nn^{}}^\lambda (𝐱)`$ $`=i\frac{1}{2}\left(J_{nn^{}}^\lambda (𝐱)[J_{n^{}n}^\lambda (𝐱)]^{}\right),`$ (78)
$`Y_{nn^{}}^{\lambda \mu }(𝐱)`$ $`=\frac{1}{2}R_{nn^{}}^{\{\lambda \mu \}}(𝐱),`$ (79)
are renormalized versions of the $`\delta `$ and $`\delta ^{}`$ mixing potentials considered previously in Eqs. (58) and (59). However, the other two terms

$`\mathrm{\Gamma }_{nn^{}}^{\lambda \mu }(𝐱)`$ $`=i\left(M_{nn^{}}^{\{\lambda \mu \}}(𝐱)[M_{n^{}n}^{\{\lambda \mu \}}(𝐱)]^{}\right),`$ (80)
$`\mathrm{\Phi }_{nn^{}}^{\lambda \mu }(𝐱)`$ $`=iR_{nn^{}}^{[\lambda \mu ]}(𝐱),`$ (81)
were not present in Eqs. (58) and (59). The term $`\mathrm{\Phi }_{nn^{}}^{\lambda \mu }(𝐱)`$, which is antisymmetric in $`\lambda `$ and $`\mu `$, is just a generalized Rashba coefficient Rashba (1960); Vas’ko (1979); Winkler (2000); Winkler et al. (2002); Winkler (2003) for multi-band Hamiltonians. However, the term $`\mathrm{\Gamma }_{nn^{}}^{\lambda \mu }(𝐱)`$, which is symmetric in $`\lambda `$ and $`\mu `$, not (g) has received little attention in the literature. Its physical interpretation will be discussed below in Sec. VIII.

Since the functions $`Z`$, $`Y`$, $`\mathrm{\Gamma }`$, and $`\mathrm{\Phi }`$ behave to lowest order as step functions at an abrupt junction, the Hamiltonian (74) shows explicitly that these functions produce interface terms proportional to $`\delta `$ or $`\delta ^{}`$.

The remaining terms ($`\widehat{\pi }`$, $`\widehat{Z}`$, etc.) in Eq. (74) that appear in front of $`\phi (𝐱)`$ are defined by the obvious generalizations of Eqs. (75)–(81). The interpretation of these terms parallels that of the terms already discussed, except that the contributions from $`\phi `$ are not as well localized at the interface.

Here it is worth noting that in the $`GW`$ approximation, Hedin (1965); Hedin and Lundqvist (1969) one has $`\mathrm{\Lambda }_{nn^{}}=\delta _{nn^{}}`$ and $`\widehat{\pi }_{nn^{}}^\lambda =\widehat{D}_{nn^{}}^{\lambda \mu }=0`$, so the potential $`\phi (𝐱)`$ does not contribute to the momentum matrix or the effective-mass tensor. In LDA, these simplifications are also valid, and one has in addition $`\widehat{Z}_{nn^{}}^\lambda =0`$, since the exchange-correlation potential is short-ranged.

For isovalent zinc-blende systems, the leading term in $`\phi ^{(1)}`$ is an octopole potential, so the contributions from the second-rank tensors $`\widehat{D}`$, $`\widehat{Y}`$, $`\widehat{\mathrm{\Gamma }}`$, and $`\widehat{\mathrm{\Phi }}`$ are negligible. These terms are non-negligible only for the linear interface dipole term in a heterovalent zinc-blende system (or for a slowly varying external potential, which is not considered explicitly here). Likewise, $`\phi ^{(2)}`$ is negligible in all terms except $`\mathrm{\Lambda }_{nn^{}}`$.

Equation (74) is the main result obtained in this paper. The qualitative form of this Hamiltonian is very similar to the Leibler Hamiltonian Leibler (1975, 1977) for slowly graded heterostructures, as amended by Takhtamirov and Volkov. Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999b) The differences are primarily due to the use of atomic pseudopotentials (rather than a model based on periodic bulk potentials), the inclusion of long-range Coulomb potentials, and the use of linear response theory to simplify the interface Hamiltonian, as discussed in Sec. IX.

The explicit form of the various material parameters in Eq. (74) for semiconductors with the zinc-blende structure is given below in Sec. VIII. First, however, the possibility of representing the material parameters as piecewise constant is considered.

## VII Simplified mathematical representation of heterostructure material properties

This section discusses several ways in which the mathematical description of material properties can be simplified. The first is to label the materials in terms of bulk compounds (e.g., GaAs) rather than atoms; the second is to approximate the properties of an abrupt junction using piecewise constant material parameters with $`\delta `$ functions and their derivatives at interfaces. For simplicity, only two-dimensional isovalent systems are considered here (see Sec. V.4).

### VII.1 Transformation from atomic to bulk-crystal description

All of the linear-response terms can be transformed immediately to a bulk-crystal representation similar to that described above in Eq. (62):

$`W_{nn^{}}^{(1)}(𝐱)`$ $`=\underset{\alpha }{{\displaystyle ^{}}}\theta ^\alpha (𝐱)W_{nn^{}}^\alpha `$ (82a)
$`=\underset{l}{{\displaystyle ^{}}}\theta ^l(𝐱)W_{nn^{}}^l,`$ (82b)
in which $`l`$ labels the different bulk materials, and \[cf. Eq. (61)\]
$$\underset{l}{}\theta ^l(𝐱)=0.$$
(82c)
For example, in a GaAs/AlAs heterostructure, the linear potentials for the two bulk media are defined by

$$\begin{array}{cc}\hfill W_{nn^{}}^{\text{GaAs}}& =W_{nn^{}}^{\text{Ga}}+W_{nn^{}}^{\text{As}},\hfill \\ \hfill W_{nn^{}}^{\text{AlAs}}& =W_{nn^{}}^{\text{Al}}+W_{nn^{}}^{\text{As}}.\hfill \end{array}$$
(83)
If the reference crystal is chosen to be GaAs, then $`W_{nn^{}}^{\text{Ga}}=W_{nn^{}}^{\text{As}}=0`$. \[This choice is made for clarity of exposition; the final results given in Eqs. (85), (90), (95), and (97) do not depend on the choice of reference crystal.\] The transformation from (82a) to (82b) is then simply

$$W_{nn^{}}^{(1)}(𝐱)=\theta ^{\text{Al}}(𝐱)W_{nn^{}}^{\text{Al}}=\theta ^{\text{Al}}(𝐱)W_{nn^{}}^{\text{AlAs}}.$$
(84)
Hence

$$\theta ^{\text{AlAs}}=\theta ^{\text{Al}},\theta ^{\text{GaAs}}=\theta ^{\text{Ga}}.$$
(85)
For a no-common-atom system such as InAs/GaSb, four different bulk potentials can be defined:

$$\begin{array}{cc}\hfill W_{nn^{}}^{\text{InAs}}& =W_{nn^{}}^{\text{In}}+W_{nn^{}}^{\text{As}},\hfill \\ \hfill W_{nn^{}}^{\text{InSb}}& =W_{nn^{}}^{\text{In}}+W_{nn^{}}^{\text{Sb}},\hfill \\ \hfill W_{nn^{}}^{\text{GaAs}}& =W_{nn^{}}^{\text{Ga}}+W_{nn^{}}^{\text{As}},\hfill \\ \hfill W_{nn^{}}^{\text{GaSb}}& =W_{nn^{}}^{\text{Ga}}+W_{nn^{}}^{\text{Sb}}.\hfill \end{array}$$
(86)
If the reference crystal is chosen to be InAs, then $`W_{nn^{}}^{\text{In}}=W_{nn^{}}^{\text{As}}=0`$. In the atomic description, the linear response potential is

$$W_{nn^{}}^{(1)}(𝐱)=\theta ^{\text{Ga}}(𝐱)W_{nn^{}}^{\text{Ga}}+\theta ^{\text{Sb}}(𝐱)W_{nn^{}}^{\text{Sb}},$$
(87)
which can be rewritten (using $`\theta ^{\text{In}}+\theta ^{\text{Ga}}=0`$ and $`\theta ^{\text{As}}+\theta ^{\text{Sb}}=0`$) in several different ways, two of which are

$$W_{nn^{}}^{(1)}=\theta ^{\text{Ga}}W_{nn^{}}^{\text{GaSb}}+(\theta ^{\text{In}}+\theta ^{\text{Sb}})W_{nn^{}}^{\text{InSb}}$$
(88a)
and
$$W_{nn^{}}^{(1)}=\theta ^{\text{Sb}}W_{nn^{}}^{\text{GaSb}}+(\theta ^{\text{Ga}}+\theta ^{\text{As}})W_{nn^{}}^{\text{GaAs}}.$$
(88b)
Equation (88a) is useful for describing a (001) heterojunction with an InSb-like interface:

$$\mathrm{}\text{–As–In–As–In–Sb–Ga–Sb–Ga–}\mathrm{},$$
(89a)
while Eq. (88b) is useful for describing a heterojunction with a GaAs-like interface:
$$\mathrm{}\text{–In–As–In–As–Ga–Sb–Ga–Sb–}\mathrm{}.$$
(89b)
In the first case one can identify the bulk functions

$`\theta ^{\text{GaSb}}`$ $`=\theta ^{\text{Ga}},`$ $`\theta ^{\text{InAs}}`$ $`=\theta ^{\text{As}},`$
$`\theta ^{\text{InSb}}`$ $`=\theta ^{\text{In}}+\theta ^{\text{Sb}},`$ $`\theta ^{\text{GaAs}}`$ $`=0,`$ (90a)
while in the second case
$`\theta ^{\text{GaSb}}`$ $`=\theta ^{\text{Sb}},`$ $`\theta ^{\text{InAs}}`$ $`=\theta ^{\text{In}},`$
$`\theta ^{\text{GaAs}}`$ $`=\theta ^{\text{Ga}}+\theta ^{\text{As}},`$ $`\theta ^{\text{InSb}}`$ $`=0.`$ (90b)
This type of transformation can be used for any term in the linear response. For the quadratic response, such a description is not appropriate, but all contributions (including the interface dipole potential) can be approximated as abrupt step functions, as discussed in Secs. V.3.1 and V.4.

### VII.2 Piecewise constant material parameters

The next step in simplifying the description of the material properties is to approximate an ideal heterostructure as piecewise constant. As a specific example, the case of a (001) heterojunction between semiconductors with the zinc-blende structure is considered here.

For this case, a convenient slab-adapted Heine (1963); Kleinman (1981) unit cell is defined by the basis vectors

$$𝐚_1=\frac{a}{2}(1,1,0),𝐚_2=\frac{a}{2}(1,1,0),𝐚_3=\frac{a}{2}(1,0,1).$$
(91)
Periodic boundary conditions are applied over the crystal volume $`\mathrm{\Omega }=𝐋_1(𝐋_2\times 𝐋_3)`$, where $`𝐋_i=N_i𝐚_i`$ and $`N_i`$ is an integer (thus $`\mathrm{\Omega }=N\mathrm{\Omega }_0`$, where $`N=N_1N_2N_3`$).

The bulk properties of the heterojunction are to be represented in terms of the periodic step function

$$\mathrm{\Theta }(𝐱)=\mathrm{\Theta }(z)=\{\begin{array}{cc}1,\hfill & 0<z<\frac{1}{2}L_z\hfill \\ 0,\hfill & \frac{1}{2}L_z<z<0\hfill \end{array}$$
(92)
in which $`L_z=\frac{1}{2}N_3a`$ is the period in the $`z`$ direction, and $`\mathrm{\Theta }(𝐱)=\mathrm{\Theta }(𝐱+𝐋_i)`$. The interface properties are to be represented by the derivatives $`\mathrm{\Theta }_z(𝐱)=\mathrm{\Theta }/z`$ and $`\mathrm{\Theta }_{zz}(𝐱)=^2\mathrm{\Theta }/z^2`$, which are periodic arrays of $`\delta `$ and $`\delta ^{}`$ functions.

The first example of an actual (001) heterojunction to be considered is a common-atom GaAs/AlAs junction. The coordinate origin is chosen to be an interface As atom, with unit-cell basis vectors $`𝝉_a=\mathrm{𝟎}`$ and $`𝝉_c=\frac{1}{4}a(1,1,1)`$ for anions and cations, respectively. The discrete function $`\theta _𝐑^\alpha `$ for $`\alpha =\text{Al}`$ is therefore

$$\theta _𝐑^{\text{Al}}=\{\begin{array}{cc}1,\hfill & 0R_3\frac{1}{2}N_31\hfill \\ 0,\hfill & \frac{1}{2}N_3R_31\hfill \end{array}$$
(93)
where the dimensionless coordinate $`R_3`$ is an integer (see Sec. V.4). From the Fourier transform of this function, one finds that for small $`k`$,

$$\theta ^{\text{Al}}(𝐤)\mathrm{\Theta }(𝐤)\frac{a^2}{96}\mathrm{\Theta }_{zz}(𝐤).$$
(94)
Upon inserting this result into Eq. (84), one finds that for $`|z|<\frac{1}{2}L_z`$, the linear band offset at a GaAs/AlAs heterojunction can be represented as

$$\begin{array}{c}W_{nn^{}}^{(1)}(𝐱)=(W_{nn^{}}^{\text{AlAs}}W_{nn^{}}^{\text{GaAs}})[\mathrm{\Theta }(z)x]\hfill \\ \hfill \frac{a^2}{96}(W_{nn^{}}^{\text{AlAs}}W_{nn^{}}^{\text{GaAs}})\delta ^{}(z),\end{array}$$
(95)
which has been written in a general form suitable for a virtual reference crystal $`\text{(GaAs)}_{1x}\text{(AlAs)}_x`$. Thus, if the band offset is treated as piecewise constant, there is an additional correction proportional to $`\delta ^{}(z)`$ that has the effect of renormalizing the matrix $`Y`$ in the Hamiltonian (74). This merely reflects the fact that $`\theta ^{\text{Al}}(z)`$ is a smooth (macroscopically averaged) step function, and the difference between a smooth step function and an abrupt step function is, to leading order, proportional to the macroscopic average of $`\delta ^{}(z)`$.

The remaining material-dependent parameters in Eq. (74) can be treated using the same approach, but for these parameters the term proportional to $`\mathrm{\Theta }_{zz}`$ yields a correction of order $`(\overline{k}a)^3\mathrm{\Delta }\overline{V}`$ or higher, which can be neglected. Thus for an ideal GaAs/AlAs junction, all of the material parameters in (74) except $`W_{nn^{}}^{(1)}(𝐱)`$ can be replaced by abrupt step functions.

For the case of a no-common-atom InAs/GaSb junction, the situation is more complicated. The example considered here is the GaAs-like junction in Eq. (89b), where the origin is the midpoint of an As–Ga bond and $`𝝉_c=𝝉_a=\frac{1}{8}a(1,1,1)`$. In this case, the small-$`k`$ behavior of the atomic distribution functions is

$$\begin{array}{cc}\hfill \theta ^{\text{Ga}}(𝐤)& \mathrm{\Theta }(𝐤)+\frac{a}{8}\mathrm{\Theta }_z(𝐤)\frac{a^2}{384}\mathrm{\Theta }_{zz}(𝐤),\hfill \\ \hfill \theta ^{\text{Sb}}(𝐤)& \mathrm{\Theta }(𝐤)\frac{a}{8}\mathrm{\Theta }_z(𝐤)\frac{a^2}{384}\mathrm{\Theta }_{zz}(𝐤).\hfill \end{array}$$
(96)
Upon inserting these results into Eq. (88b), one finds that for $`|z|<\frac{1}{2}L_z`$, the linear band offset may be approximated by

$$\begin{array}{c}W_{nn^{}}^{(1)}(𝐱)=(W_{nn^{}}^{\text{GaSb}}W_{nn^{}}^{\text{InAs}})[\mathrm{\Theta }(z)x]\hfill \\ \hfill +\frac{a}{8}(2W_{nn^{}}^{\text{GaAs}}W_{nn^{}}^{\text{InAs}}W_{nn^{}}^{\text{GaSb}})\delta (z)\\ \hfill \frac{a^2}{384}(W_{nn^{}}^{\text{GaSb}}W_{nn^{}}^{\text{InAs}})\delta ^{}(z),\end{array}$$
(97a)
which is written in a form suitable for the virtual reference crystal $`\text{(InAs)}_{1x}\text{(GaSb)}_x`$. Note that this has the form
$$\begin{array}{c}W_{nn^{}}^{(1)}(𝐱)=(W_{nn^{}}^{(+)}W_{nn^{}}^{()})[\mathrm{\Theta }(z)x]\hfill \\ \hfill +\frac{a}{8}(2W_{nn^{}}^{\text{(i)}}W_{nn^{}}^{(+)}W_{nn^{}}^{()})\delta (z)\\ \hfill \frac{a^2}{384}(W_{nn^{}}^{(+)}W_{nn^{}}^{()})\delta ^{}(z),\end{array}$$
(97b)
in which $`(+)`$ and $`()`$ label the bulk materials to the right and left of the junction, (i) labels the “interface” material, and the reference crystal is $`(+)_{1x}()_x`$. This result is valid for any no-common-atom (001) junction, regardless of the convention chosen for the $`𝝉`$ vectors. Thus for a no-common-atom junction, $`W_{nn^{}}^{(1)}(𝐱)`$ contains terms that renormalize both the $`Z`$ and $`Y`$ matrices.

The $`\delta (z)`$ term in (97) has a very simple physical interpretation in which half a monolayer ($`a/4`$) at the interface is occupied by GaAs instead of either InAs or GaSb. This concept has been used previously in the construction of envelope-function models for electrons Waterman et al. (1993) and phonons. Foreman (1995b)

Note that if one replaces Sb with As in Eq. (97a), the result does not reduce directly to an expression of the form (95). However, the difference is merely due to the different choice of coordinate origin in the two cases. To leading order, $`\mathrm{\Theta }(z\frac{1}{8}a)\mathrm{\Theta }(z)\frac{1}{8}a\delta (z)+\frac{1}{128}a^2\delta ^{}(z)`$, $`\delta (z\frac{1}{8}a)\delta (z)\frac{1}{8}a\delta ^{}(z)`$, and $`\delta ^{}(z\frac{1}{8}a)\delta ^{}(z)`$. Thus, replacing Sb with As and shifting the origin by $`\frac{1}{8}a(1,1,1)`$ does indeed reduce (97a) to (95).

The result (97b) can be applied to all other linear-response terms in the Hamiltonian (74), but since the $`\mathrm{\Delta }\pi (𝐱)`$ and $`Z(𝐱)`$ terms are already of order $`(\overline{k}a)\mathrm{\Delta }\overline{V}`$, the $`\delta ^{}(z)`$ term in (97b) yields a negligible correction of order $`(\overline{k}a)^3\mathrm{\Delta }\overline{V}`$. Thus

$$\begin{array}{c}Z_{nn^{}}^\lambda (𝐱)=(Z_{nn^{}}^{\lambda (+)}Z_{nn^{}}^{\lambda ()})[\mathrm{\Theta }(z)x]\hfill \\ \hfill +\frac{a}{8}(2Z_{nn^{}}^{\lambda \text{(i)}}Z_{nn^{}}^{\lambda (+)}Z_{nn^{}}^{\lambda ()})\delta (z),\end{array}$$
(98)
and we see that $`Z`$ generates a correction to $`Y`$, while $`\pi `$ generates a correction to $`\mathrm{\Gamma }`$. However, the position dependence of all remaining terms in (74) (including the quadratic response) can be represented as a simple step function:

$$Y_{nn^{}}^{\lambda \mu }(𝐱)=(Y_{nn^{}}^{\lambda \mu (+)}Y_{nn^{}}^{\lambda \mu ()})[\mathrm{\Theta }(z)x],$$
(99)
since these terms are already of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$.

## VIII Theory of invariants

In this section the method of invariants Luttinger (1956); Suzuki and Hensel (1974); Bir and Pikus (1974); Trebin et al. (1979); Ivchenko and Pikus (1997) is used to construct the explicit form of the Hamiltonian (74) for $`\mathrm{\Gamma }`$ electrons in semiconductors with the zinc-blende structure. Parmenter (1955); Dresselhaus (1955) The results are then compared with the interface Hamiltonian given by the method of invariants for a common-atom (001) junction, including all (linear and nonlinear) interface terms of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$ that are permitted by symmetry. Finally, the interface terms of order $`(\overline{k}a)\mathrm{\Delta }\overline{V}`$ are considered for the case of $`\mathrm{\Gamma }`$$`X`$ coupling in GaAs/AlAs. The results in this section make use of the time reversal and crystal symmetry properties of the self-energy that were derived in Ref. Foreman, 2005b.

### VIII.1 Generalized Leibler Hamiltonian

As demonstrated above, to within terms of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$, the position dependence of the interface parameters can be calculated entirely in terms of bulk-like matrix elements modulated by the atomic form factors $`\theta ^\alpha (𝐱)`$. Thus, the Hamiltonian for a zinc-blende heterostructure of arbitrary composition can be determined by constructing invariants transforming as $`\mathrm{\Gamma }_1`$ under the symmetry operations of the $`T_d`$ group. This amounts to treating $`\theta ^\alpha (𝐱)`$ as a “slowly varying” function that is an invariant of $`T_d`$. Leibler (1977)

The relevant basis functions for the representations of $`T_d`$ are given in Table 1.

The specific example to be considered here is that of the $`\mathrm{\Gamma }_8`$ valence band, with $`\mathrm{\Gamma }_{15}`$ and $`\mathrm{\Gamma }_6`$ derivable as special cases of the $`\mathrm{\Gamma }_8`$ results. The extension of these results to multi-band Hamiltonians (e.g., $`\mathrm{\Gamma }_6\mathrm{\Gamma }_7\mathrm{\Gamma }_8`$) can be handled using the methods of Refs. Suzuki and Hensel, 1974; Bir and Pikus, 1974; Trebin et al., 1979, but is not considered explicitly here.

The momentum matrix $`\pi ^\lambda (𝐱)`$ must transform as a vector ($`\mathrm{\Gamma }_{15}`$) that is odd under time reversal, which leaves only one possibility:

$$\pi ^\lambda (𝐱)=V_\lambda \xi (𝐱),$$
(100)
in which $`V_x\{J_x(J_y^2J_z^2)\}`$ is a product of angular momentum matrices $`J_\nu `$ for a particle with spin $`\frac{3}{2}`$, and $`\xi (𝐱)`$ is some linear combination of $`\theta ^l(𝐱)`$ functions. The corresponding results for $`\mathrm{\Gamma }_{15}`$ and $`\mathrm{\Gamma }_6`$ states are obtained by replacing $`J_\nu I_\nu `$ (spin 1) and $`J_\nu \sigma _\nu `$ (spin $`\frac{1}{2}`$), respectively. This yields zero in both cases, hence $`\pi ^\lambda (𝐱)`$ occurs only for $`\mathrm{\Gamma }_8`$ and is a relativistic effect.

For the effective-mass tensor $`D^{\lambda \mu }(𝐱)`$ one obtains the well-known result Luttinger (1956)

$$\begin{array}{c}D^{\lambda \mu }(𝐱)=\frac{1}{2}\gamma _1(𝐱)\delta _{\lambda \mu }1\gamma _2(𝐱)\delta _{\lambda \mu }(J_\lambda ^2\frac{1}{3}J^2)\hfill \\ \hfill \gamma _3(𝐱)(1\delta _{\lambda \mu })\{J_\lambda J_\mu \}\\ \hfill +iϵ_{\lambda \mu \nu }[\kappa (𝐱)J_\nu +q(𝐱)J_\nu ^3],\end{array}$$
(101)
in which the Luttinger parameters $`\gamma _1(𝐱)`$, $`\gamma _2(𝐱)`$, $`\gamma _3(𝐱)`$, $`\kappa (𝐱)`$, and $`q(𝐱)`$ are all linear combinations of $`\theta ^l(𝐱)`$ functions, and $`ϵ_{\lambda \mu \nu }`$ is the antisymmetric unit tensor. For $`\mathrm{\Gamma }_{15}`$ the parameter $`q`$ is not independent (since $`I_\nu ^3=I_\nu `$), whereas for $`\mathrm{\Gamma }_6`$ only $`\gamma _1`$ and $`\kappa `$ are independent.

The $`Z`$ matrix has the same symmetry as the coordinate operator (i.e., a vector that is even under time reversal), so it has the form Ivchenko et al. (1996)

$$Z^\lambda (𝐱)=|ϵ_{\lambda \mu \nu }|\{J_\mu J_\nu \}\zeta (𝐱),$$
(102)
where $`\zeta (𝐱)`$ is a linear combination of $`\theta ^l(𝐱)`$ functions. This result has the same form for $`\mathrm{\Gamma }_{15}`$, but vanishes for $`\mathrm{\Gamma }_6`$. The coupling (102) generates a zone-center mixing of heavy and light holes, and was proposed independently in Refs. Aleĭner and Ivchenko, 1992, Foreman, 1995a, Foreman, 1996, and Krebs and Voisin, 1996.

The $`Y`$ interface matrix has the same symmetry as the symmetric part of $`D`$, so it can be written

$$\begin{array}{c}Y^{\lambda \mu }(𝐱)=\frac{1}{2}\eta _1(𝐱)\delta _{\lambda \mu }1\eta _2(𝐱)\delta _{\lambda \mu }(J_\lambda ^2\frac{1}{3}J^2)\hfill \\ \hfill \eta _3(𝐱)(1\delta _{\lambda \mu })\{J_\lambda J_\mu \}.\end{array}$$
(103)
This term has not been studied previously for $`\mathrm{\Gamma }_8`$ states, although the corresponding term for $`\mathrm{\Gamma }_6`$ electrons (a direct analog of the Darwin term from the Dirac equation) is well known. Leibler (1975, 1977); Young (1989)

The $`\mathrm{\Gamma }`$ matrix has the form

$$\mathrm{\Gamma }^{\lambda \mu }(𝐱)=|ϵ_{\lambda \mu \nu }|V_\nu \chi (𝐱),$$
(104)
in which $`\chi (𝐱)`$ is a linear combination of $`\theta ^l(𝐱)`$ functions. The closely related coupling $`\widehat{\mathrm{\Gamma }}^{\lambda \mu }`$ arising from an external electric field was proposed recently in Ref. Ganichev et al., 2002b as a possible mechanism for the linear and circular photogalvanic effects. This term is relativistic in origin and does not occur for $`\mathrm{\Gamma }_{15}`$ or $`\mathrm{\Gamma }_6`$.

The matrix $`\mathrm{\Phi }`$ has the same symmetry as the antisymmetric part of $`D`$, except that it is odd rather than even under time reversal. It can therefore be written as

$$\mathrm{\Phi }^{\lambda \mu }(𝐱)=ϵ_{\lambda \mu \nu }[\kappa ^{}(𝐱)J_\nu +q^{}(𝐱)J_\nu ^3].$$
(105)
For $`\mathrm{\Gamma }_6`$ electrons (where $`q^{}`$ is not independent) this is the analog of the Rashba spin-splitting effect for low-symmetry bulk semiconductors, Rashba (1960) which occurs also in heterostructures of cubic semiconductors Leibler (1977); Vas’ko (1979); Bastard (1981, 1982) due to the reduced symmetry at a surface or interface. The corresponding Hamiltonian for $`\mathrm{\Gamma }_8`$ states was proposed in Ref. Foreman, 1993, and has received renewed attention as the valence-band Rashba coupling Winkler (2000); Winkler et al. (2002); Winkler (2003); Bernevig and Zhang (2005) in recent years.

The coefficients of the various terms involving $`\phi `$ in Eq. (74) have the same form as those already given. The symmetry restrictions on the Taylor series expansions for $`n^{(1)}`$ and $`n^{(2)}`$ were given in Ref. Foreman, 2005b.

### VIII.2 GaAs/AlAs (001) heterojunction

The significance of these results is now investigated by comparing them with all interface terms of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$ allowed by symmetry for a GaAs/AlAs (001) heterojunction, which has the point group $`C_{2v}`$. To better understand this comparison, it is helpful to begin by studying the symmetry properties of the linear and quadratic response in this system.

The starting point is the observation that although the point group of the heterojunction is $`C_{2v}`$ (with the coordinate origin at an interface As atom), if the reference crystal is chosen to be the virtual crystal Al0.5Ga0.5As, then the perturbation due to the ionic pseudopotentials has a higher symmetry: it transforms according to the representation $`X_3`$ of the $`D_{2d}`$ group. In other words, it transforms according to the identity representation $`\mathrm{\Delta }_1`$ of the $`C_{2v}`$ group (none of whose operations change the $`z`$ coordinate), but in addition it has odd parity with respect to those elements of $`D_{2d}`$ that change the sign of $`z`$. This occurs because the atomic mole fraction $`\theta _𝐑^\alpha `$ itself transforms as $`X_3(D_{2d})`$, whereas the atomic perturbation $`\mathrm{\Delta }v_{\text{ion}}^\alpha (𝐱,𝐱^{})`$ has the site symmetry $`\mathrm{\Gamma }_1(T_d)`$ or $`X_1(D_{2d})`$. Therefore, $`\mathrm{\Delta }V_{\text{ion}}(𝐱,𝐱^{})`$ transforms as $`X_1X_3=X_3`$ under the operations of $`D_{2d}`$.

To determine the behavior of the screened potential $`\mathrm{\Delta }V(𝐱,𝐱^{})`$, note that since the vertex functions Foreman (2005b) $`\mathrm{\Gamma }^{(1)}`$ and $`\mathrm{\Gamma }^{(2)}`$ of the reference crystal transform as $`\mathrm{\Gamma }_1(T_d)`$, the linear response $`V^{(1)}(𝐱,𝐱^{})`$ transforms as $`X_1X_3=X_3`$, while the quadratic response $`V^{(2)}(𝐱,𝐱^{})`$ transforms as $`X_1X_3X_3=X_1`$. (In general, all odd-order terms in the response transform as $`X_3`$, while all even-order terms transform as $`X_1`$.) Thus, the total response $`\mathrm{\Delta }V(𝐱,𝐱^{})=V^{(1)}(𝐱,𝐱^{})+V^{(2)}(𝐱,𝐱^{})`$ transforms as neither $`X_3`$ nor $`X_1`$, but as $`\mathrm{\Delta }_1(C_{2v})`$.

If the quadratic response is of the same order as the linear response, then one must include all possible invariants of $`C_{2v}`$ that are of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$ when constructing the interface Hamiltonian. Basis functions for $`D_{2d}`$ and $`C_{2v}`$ are given in Tables 2 and 3, respectively; note that both $`X_1`$ and $`X_3`$ are compatible with $`\mathrm{\Delta }_1`$.

All interface invariants of $`C_{2v}`$ that are (1) hermitian, (2) time-reversal invariant, and (3) of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$ or less are listed in Table 4.

The second column in this table lists the symmetry of each term under the operations of $`D_{2d}`$. The third column indicates whether the term can be constructed as an invariant of $`T_d`$, treating $`\mathrm{\Theta }(z)`$ as an invariant of $`T_d`$. (As discussed above, the latter approach can be used for all linear-response interface terms.) Note that in the $`T_d`$ symmetry analysis, $`\delta (z)=d\mathrm{\Theta }/dz`$ transforms as $`z`$ ($`\mathrm{\Gamma }_{15}`$) and $`\delta ^{}(z)`$ transforms as $`z^2`$ ($`\mathrm{\Gamma }_1`$ or $`\mathrm{\Gamma }_{12}`$), whereas in the $`D_{2d}`$ symmetry analysis, $`\delta (z)`$ transforms as $`z^2`$ ($`X_1`$) and $`\delta ^{}(z)`$ transforms as $`z`$ ($`X_3`$).

In agreement with the general symmetry properties of $`V^{(1)}`$ and $`V^{(2)}`$ derived earlier, Table 4 shows that all linear interface terms \[derivable from $`\mathrm{\Gamma }_1(T_d)`$\] transform as $`X_3(D_{2d})`$. The remaining terms that transform as $`X_1`$ all originate in the quadratic or higher-order response. Thus, if one accounts for the smallness of the quadratic response (Sec. III.2), seven out of the thirteen possible interface invariants for (001) GaAs/AlAs heterojunctions can be omitted because they are actually of order $`(\overline{k}a)^3\mathrm{\Delta }\overline{V}`$ or higher. This represents a major simplification over the general case in which all response terms are of the same order.

For a no-common-atom heterojunction, the results of Sec. VIII.1 remain valid (as they depend only on the bulk symmetry), but the linear response now has only the symmetry $`\mathrm{\Delta }_1(C_{2v})`$. Thus the terms labeled $`X_1`$ in Table 4 now arise (in general) even in linear-response theory, due to the renormalization effects described above in Eqs. (97) and (98). For example, the $`\delta (z)`$ term in Eq. (97) generates the term in the first row of Table 4. However, the term in the second row of Table 4 is still zero (within linear response), because $`W_{nn^{}}^{(1)}(𝐱)`$ in Eq. (97) couples only states of the same symmetry, with $`W_{nn^{}}^{(1)}(𝐱)=\delta _{nn^{}}W_{nn}^{(1)}(𝐱)`$ and $`W_{nn}^{(1)}(𝐱)=W_{n^{}n^{}}^{(1)}(𝐱)`$ whenever $`E_n=E_n^{}`$.

### VIII.3 $`\mathrm{\Gamma }`$$`X`$ coupling

If one extends the above analysis to the case of intervalley $`\mathrm{\Gamma }`$$`X`$ coupling at a (001) GaAs/AlAs junction, it is immediately apparent (because the $`Z`$ matrix has the symmetry of a coordinate matrix) that the linear response produces $`\delta (z)`$ coupling between the bands

$$\mathrm{\Gamma }_1\text{}X_{3z},X_{1z}\text{}X_{3z},X_{3x}\text{}X_{3y},$$
(106)
but no $`\delta (z)`$ coupling between the bands

$$\mathrm{\Gamma }_1\text{}X_{1z},X_{1x}\text{}X_{1y}.$$
(107)
Coupling of the latter type only occurs in the quadratic response, or in the linear response from terms proportional to $`\delta ^{}(z)`$ or $`i\{P_z\delta (z)\}`$. These conclusions agree with those of Ref. Foreman, 1998, which were derived from a model potential constructed from a linear superposition of atomic-like pseudopotentials.

Note, however, that these results hold only for an ideal heterojunction. If interdiffusion of Ga and Al atoms breaks the $`X_3`$ symmetry of the linear response, then there will be $`\delta (z)`$ coupling between the bands (107) due to terms similar to those derived in Eq. (97).

## IX Discussion

Since a summary of the main results of this paper has already been given in Sec. I.3, this section will be limited to a discussion of the differences between the present theory and previous envelope-function theories that have appeared in the literature. The principle differences are that (1) the present Hamiltonian is constructed from atomic pseudopotentials rather than the periodic potential of a bulk crystal; (2) the present theory is self-consistent, accounting fully for electron-electron Coulomb interactions; and (3) the present approach uses linear and quadratic response theory to simplify the functional form of the heterostructure Hamiltonian.

Starting with point (1), previous envelope-function theories Leibler (1975, 1977); Burt (1988, 1992, 1994, 1999); Young (1989); Cuypers and van Haeringen (1991); Karavaev and Tikhodeev (1991); Karavaev and Krivorotov (1996a, b); Foreman (1993, 1995a, 1996, 1998, 2001); Burt and Foreman (1999); Klipstein (1999); Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b); Cortez et al. (2000, 2001); Smith and Mailhiot (1986, 1990); Mailhiot and Smith (1990); Trzeciakowski (1988); Ando et al. (1989); Cuypers and van Haeringen (1992); Cuypers and van Haeringen (1993a, b); Aversa and Sipe (1993, 1994); Grinyaev et al. (1994) were based on model potentials of the form (62), which is a linear combination of periodic bulk potentials multiplied by step-like functions $`\theta ^l(𝐱)`$. The form of $`\theta ^l(𝐱)`$ near an interface was either specified as part of the model (e.g., abrupt step function) or treated as unknown—and in principle unknowable, at least within the model potential approach.

Although the present theory led to a similar expression (82b), this was a derived result based on linear response theory, in which the functions $`\theta ^l(𝐱)`$ have a known form determined by the distribution $`\theta _𝐑^\alpha `$ of atoms in the heterostructure. The potential (82b) includes contributions from the “interface” materials (e.g., GaAs and InSb in an InAs/GaSb heterostructure), which are often omitted in the model-potential approach. Furthermore, for the model potential, the coefficients of the $`\delta `$ terms in the Hamiltonian depend on (i) the values of the bulk potentials at the reciprocal lattice vectors $`𝐆`$; and (ii) whether the functions $`\theta ^l(𝐱)`$ are smooth (zero outside the Brillouin zone) or sharp, with some $`\delta `$ terms vanishing in the smooth case. Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b) In the present theory these coefficients are determined not by $`\theta ^l(𝐱)`$, but instead by the properties of the linear atomic pseudopotentials in a finite neighborhood of each $`𝐆`$.

For point (2), the present theory includes long-range multipole Coulomb potentials (because the self-energy is nonanalytic at $`𝐪=𝐆`$), which were not considered in previous studies. These terms have no qualitative effect in two-dimensional systems (other than to renormalize the band offsets), but they are present in quantum wires and dots, even for isovalent heterostructures.

Finally, for point (3), if the functions $`\theta ^l(𝐱)`$ are treated as unknown in the model potential approach, one has no way of knowing how large the various interface terms are, and all terms permitted by symmetry should in principle be included. Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b) However, since the linear response is dominant (for typical heterostructures) in a first-principles theory, one can eliminate all nonlinear interface band-mixing terms, thereby simplifying the Hamiltonian considerably.

Aside from these differences, the qualitative form of the Hamiltonian derived here is very similar to those derived by Leibler Leibler (1975, 1977) for slowly graded heterostructures and Takhtamirov and Volkov Takhtamirov and Volkov (1997); Volkov and Takhtamirov (1997); Takhtamirov and Volkov (1999a, b, 2000b) for abrupt heterostructures. The basic structure of these Hamiltonians is similar because the same type of perturbation theory Luttinger and Kohn (1955); Bir and Pikus (1974); Winkler (2003) was used in their derivation.

Other differences arise in theories using different approximation techniques. The well-known theory of Burt Burt (1992) is based on the Luttinger-Kohn representation Luttinger and Kohn (1955); Leibler (1975) with an energy-dependent approximation similar to Löwdin perturbation theory Löwdin (1951) used to eliminate the coupling to remote bands. Estimation of the error involved in this approximation shows that Burt’s general theory Burt (1992) includes all terms of order $`(\overline{k}a)\mathrm{\Delta }\overline{V}`$, but retains some terms of order $`(\overline{k}a)^2\mathrm{\Delta }\overline{V}`$ while omitting others. In particular, his results give an essentially correct (within the model potential approach) description of the $`\delta `$-function term $`Z`$ and the Rashba term $`\mathrm{\Phi }`$, but omit the contributions from remote-band coupling to the terms $`Y`$ and $`\mathrm{\Gamma }`$.

In Ref. Foreman, 1993, I used Burt’s theory, but introduced an additional approximation (also adopted in Refs. Burt, 1999 and Foreman, 2001) whereby the energy eigenvalue in the denominator of the energy-dependent effective mass was replaced with the position-dependent energy of the bulk valence band maximum. This has the effect of replacing the Rashba coefficient in Eq. (81) with $`iD_{nn^{}}^{[\lambda \mu ]}(𝐱)`$. In other words, the Rashba parameters $`\kappa ^{}`$ and $`q^{}`$ of Eq. (105) were replaced by the Luttinger parameters $`\kappa `$ and $`q`$ of Eq. (101). This simplifies applications of the theory because $`\kappa `$ and $`q`$ are known from bulk magnetoabsorption measurements, Pidgeon and Brown (1966); Pidgeon and Groves (1969) and $`\kappa `$ can also be estimated from $`\gamma _1`$, $`\gamma _2`$, and $`\gamma _3`$. Pidgeon and Brown (1966); Foreman (1993) However, since the true Rashba parameters $`\kappa ^{}`$ and $`q^{}`$ cannot be determined from bulk measurements, this simplification will produce values that have the correct order of magnitude but are quantitatively incorrect.

A similar approximation was used in Refs. Krebs and Voisin, 1996, Cortez et al., 2000, and Cortez et al., 2001 to estimate the $`\delta `$-function mixing parameter for light and heavy holes \[see Table 4 and Eq. (102)\] from the valence-band offset. The original proposal was based on the $`H_{BF}`$ model, Krebs and Voisin (1996) in which the zone-center $`\mathrm{\Gamma }_{15}`$ Hamiltonian is written as a sum of operators $`B=\frac{1}{2}+\{I_xI_y\}`$ and $`F=\frac{1}{2}\{I_xI_y\}`$ multiplied by position-dependent bulk valence-band energies $`E_v`$. For a (001) GaAs/AlAs heterojunction this yields the $`\mathrm{\Gamma }_{15}`$ interface Hamiltonian

$$\frac{a}{4}\{I_xI_y\}(E_v^{\text{GaAs}}E_v^{\text{AlAs}})\delta (z),$$
(108)
whereas for an InAs/GaSb junction with a GaAs-like interface one obtains

$$\frac{a}{4}(12\{I_xI_y\})[E_v^{\text{GaAs}}\frac{1}{2}(E_v^{\text{InAs}}+E_v^{\text{GaSb}})]\delta (z).$$
(109)
The corresponding results for $`\mathrm{\Gamma }_8`$ are given by the substitution $`\{I_xI_y\}\frac{1}{3}\{J_xJ_y\}`$. Luttinger (1956) Now the diagonal term \[proportional to $`1\delta (z)`$\] in Eq. (109) is identical (within linear response) to the $`\delta (z)`$ term derived above in Eq. (97). However, the valence-band mixing term derived in Eqs. (78) and (102) is not related to the valence-band offset. Therefore, although the model of Refs. Krebs and Voisin, 1996, Cortez et al., 2000, and Cortez et al., 2001 yields a coupling of the correct symmetry and order of magnitude, its numerical value is not reliable. Foreman (1999) The same conclusion was reached in Ref. Magri and Zunger, 2000.

In Ref. Foreman, 1996, I proposed a smooth model potential not (h) for abrupt heterostructures that yields a valence-band mixing determined chiefly by the change in Bloch functions at an interface. This result is of course not valid for general model potentials Foreman (1996, 1998) (despite claims to the contraryKaravaev and Krivorotov (1996a, b)), and it is not supported by the present work. However, in Ref. Takhtamirov and Volkov, 1999b it was argued that the valence-band mixing generated by Bloch-function differences is negligible in comparison to that generated by a discontinuous model potential. Although the latter does not appear in the present first-principles theory, its analog is the $`Z`$ term shown in Eq. (59). The former contribution is given by the term proportional to $`W_{in^{}}^{(1)}(𝐪)`$ in Eqs. (118) and (78). Inspection of these terms shows that both are of order $`(\overline{k}a)\mathrm{\Delta }\overline{V}`$. Therefore, neither the model of Ref. Foreman, 1996 \[which omits the interface contributions from Eq. (58)\] nor the model of Refs. Volkov and Takhtamirov, 1997; Takhtamirov and Volkov, 1999a, b (which omits the contribution from Bloch-function differences) gives a correct value for the valence-band mixing coefficient.

In Refs. Im et al., 1999, 2000; Klipstein, 2001, Klipstein, 1999, and Takhtamirov and Volkov, 2000b it was proposed that intervalley $`\mathrm{\Gamma }_1`$$`X_{1z}`$ and $`X_{1x}`$$`X_{1y}`$ coupling should be proportional to $`\delta (z)`$ and of the same order of magnitude as $`\mathrm{\Gamma }_1`$$`X_{3z}`$, $`X_{1z}`$$`X_{3z}`$, and $`X_{3x}`$$`X_{3y}`$ coupling. However, as shown above in Sec. VIII.3, the former types of coupling should generally be substantially weaker than the latter for ideal interfaces. The same conclusion was reached in Ref. Wang and Zunger, 1997, where it was shown that despite the very small magnitude of the direct $`\mathrm{\Gamma }_1`$$`X_{1z}`$ coupling, there is a net effective $`\mathrm{\Gamma }_1`$$`X_{1z}`$ coupling (as revealed by the anticrossing of superlattice subbands in Fig. 1 of Ref. Wang and Zunger, 1997) generated by $`\mathrm{\Gamma }_1`$$`X_{3z}`$, $`X_{1z}`$$`X_{3z}`$, and $`𝐤𝐩`$ couplings that is only a factor of two or three smaller than the net $`\mathrm{\Gamma }_1`$$`X_{3z}`$ coupling.

However, in Refs. Im et al., 1999, 2000; Klipstein, 2001 it was found that $`X_x`$$`X_y`$ mixing experiments could not adequately be explained without including $`\delta `$-function $`X_{1x}`$$`X_{1y}`$ coupling of comparable magnitude to $`X_{3x}`$$`X_{3y}`$ coupling. Nevertheless, since Refs. Im et al., 1999, 2000; Klipstein, 2001 did not test the effect of $`\delta ^{}(z)`$ $`X_{1x}`$$`X_{1y}`$ coupling, it is not yet clear that their experiments require such a large $`\delta (z)`$ coupling. If $`\delta ^{}(z)`$ or $`i\{P_z\delta (z)\}`$ coupling is unable to explain their results, then their data may indicate the presence of nonideal interfaces, as discussed in Sec. VIII.3. Another possible explanation would be spin-orbit coupling, which generates a direct $`\delta (z)`$ coupling between $`X_{6x}`$ and $`X_{6y}`$ even for ideal interfaces.

Two important physical effects not considered in this paper are the contributions from alloying (beyond the virtual crystal approximation) and from strain (in the bulk or at an interface). The former effect is in principle already encompassed by the present formalism (for lattice-matched systems in which linear response is dominant), although for practical applications one would need to perform a detailed analysis of the significance of different intervalley mixing effects. Wang and Zunger (1997) The contribution from strain would require a nontrivial extension of this theory, but since an ab initio linear-response approach to strain has already been successfully applied to In0.53Ga0.47As/InP, Baroni et al. (1989); Peressi et al. (1990) Si/Ge, Colombo et al. (1991); Peressi and Baroni (1994) and InAs/GaSb Montanari et al. (1996) within LDA, the general perturbation scheme used in this paper should work for lattice-mismatched systems also.

###### Acknowledgements.

I am grateful to Eduard Takhtamirov and Roland Winkler for valuable comments on an earlier version of the manuscript. This work was supported by Hong Kong RGC Grant HKUST6139/00P and UGC Grant HIA03/04.SC02.

## Appendix A Imaginary part of the self-energy

Luttinger Luttinger (1961) has shown for homogeneous systems that $`\mathrm{\Sigma }^{\text{(i)}}`$ vanishes at the Fermi surface to all orders in perturbation theory. The generalization of his result to the case of inhomogeneous systems is $`\mathrm{\Sigma }^{\text{(i)}}(\omega \pm i0^+)=J(\omega )`$, where $`\omega `$ is real and Negele and Orland (1998); Farid (1999)

$$J(\omega )=\{\begin{array}{cc}C_+(\omega \mu _+)^2,\hfill & \omega \mu _+\hfill \\ 0,\hfill & \mu _{}\omega \mu _+\hfill \\ C_{}(\omega \mu _{})^2,\hfill & \omega \mu _{}\hfill \end{array}$$
(110)
for $`\omega `$ near $`\mu _\pm \mu \pm \frac{1}{2}E_\text{g}`$. Here $`\mu `$ is the chemical potential (in the limit of zero temperature), $`E_\text{g}`$ is the energy gap, and $`C_\pm 0`$ is independent of $`\omega `$. The magnitude of $`C_\pm `$ may be estimated from the calculations of Quinn and Ferrell Quinn and Ferrell (1958) for a degenerate homogeneous electron gas. Their results may be expressed as $`C_\pm =1/E_\text{i}`$, where Quinn and Ferrell (1958); Callaway (1991)

$$E_\text{i}=\frac{256}{\sqrt{3}\pi ^2}\frac{ϵ_\text{F}^2}{\omega _\text{p}},$$
(111)
in which $`ϵ_\text{F}`$ is the Fermi energy and $`\omega _\text{p}`$ is the plasma frequency. If one treats semiconductors with the diamond or zinc-blende structure as a homogeneous gas with a density of eight valence electrons per primitive unit cell, one obtains values of $`E_\text{i}126`$ eV for GaAs and $`E_\text{i}140`$ eV for Si. This is very large in comparison to the energy gap and suggests that $`\mathrm{\Sigma }^{\text{(i)}}`$ will have a negligible influence on the band structure in typical semiconductor heterostructures (although it has an important qualitative effect in producing a finite quasiparticle lifetime).

Since a real semiconductor is not a degenerate homogeneous electron gas, one may question whether this estimate is reliable. However, the calculations of Fleszar and Hanke Fleszar and Hanke (1997) for Si yield $`E_\text{i}100`$ eV (estimated from Fig. 3 of Ref. Fleszar and Hanke, 1997), with some asymmetry between electrons and holes. Thus it is reasonable in heterostructures to assume that $`E_\text{i}\overline{E}_\text{g}`$. Since $`(\omega \mu _\pm )^2`$ is of order $`(\mathrm{\Delta }\overline{V})^2`$, this means that $`\mathrm{\Sigma }^{\text{(i)}}`$ is negligible under the present perturbation scheme. Thus it will, for the most part, not be considered explicitly in the envelope-function equations derived here. However, the leading contribution from $`\mathrm{\Sigma }^{\text{(i)}}`$ is noted in Eq. (56) so that it may be included if desired.

## Appendix B Cutoff function

This appendix considers possible definitions of the Brillouin zone cutoff function $`B(𝐤)`$ introduced in Eq. (60). The most obvious choice \[from Eq. (27a)\] would be

$$B(𝐤)=\{\begin{array}{cc}1,\hfill & 𝐤\mathrm{\Omega }_0^{}\hfill \\ 0,\hfill & 𝐤\mathrm{\Omega }_0^{}\hfill \end{array}$$
(112)
which was used previously in Refs. Burt, 1988 and Burt, 1992. However, the sharp cutoff at the boundary of $`\mathrm{\Omega }_0^{}`$ has the undesirable effect of producing Gibbs oscillations in $`𝐱`$ space. Other possibilities that eliminate this problem are $`B(𝐤)=[1+\beta (ka)^2]\mathrm{exp}[\beta (ka)^2]`$ and $`B(𝐤)=\mathrm{exp}[\beta (ka)^4]`$, in which $`\beta `$ is some number of order 1. These choices of $`B(𝐤)`$ are smooth, spherically symmetric, and introduce a negligible error of order $`(\overline{k}a)^4\mathrm{\Delta }\overline{V}`$.

For some applications one need not introduce any cutoff function at all, as an alternative power-series approximation for $`\theta ^\alpha (𝐤)`$ is often more convenient. This approximation is discussed in Sec. VII.2.

## Appendix C Renormalized effective-mass parameters

This appendix presents some details of the perturbation theory used in Sec. VI. For states $`m,m^{}A`$, the effective Hamiltonian $`\overline{H}`$ is given by Bir and Pikus (1974); Winkler (2003)

$$\begin{array}{c}\overline{H}_{mm^{}}=H_{mm^{}}+\frac{1}{2}\underset{i}{\overset{B}{}}H_{mi}^{}H_{im^{}}^{}\left(\frac{1}{E_{mi}}+\frac{1}{E_{m^{}i}}\right)+\frac{1}{2}\underset{i}{\overset{B}{}}\underset{j}{\overset{B}{}}H_{mi}^{}H_{ij}^{}H_{jm^{}}^{}\left(\frac{1}{E_{mi}E_{mj}}+\frac{1}{E_{m^{}i}E_{m^{}j}}\right)\hfill \\ \hfill \frac{1}{2}\underset{i}{\overset{B}{}}\underset{m^{\prime \prime }}{\overset{A}{}}\left(\frac{H_{mm^{\prime \prime }}^{}H_{m^{\prime \prime }i}^{}H_{im^{}}^{}}{E_{mi}E_{m^{\prime \prime }i}}+\frac{H_{mi}^{}H_{im^{\prime \prime }}^{}H_{m^{\prime \prime }m^{}}^{}}{E_{m^{}i}E_{m^{\prime \prime }i}}\right)+O\left((H^{})^4\right),\end{array}$$
(113)
in which $`E_{mi}=E_mE_i`$. This expression can be used to derive most of the terms in the effective-mass Hamiltonian (73). The coefficients of the quadratic, cubic, and quartic dispersion terms of the reference crystal are given by

$$D_{nn^{}}^{\lambda \mu }=\stackrel{~}{D}_{nn^{}}^{\lambda \mu }+\frac{1}{2}\underset{i}{\overset{B}{}}\pi _{ni}^\lambda \pi _{in^{}}^\mu \left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right),$$
(114)
$$\begin{array}{c}C_{nn^{}}^{\lambda \mu \kappa }=\stackrel{~}{C}_{nn^{}}^{\lambda \mu \kappa }+\frac{1}{2}\underset{i}{\overset{B}{}}(\pi _{ni}^\lambda \stackrel{~}{D}_{in^{}}^{\mu \kappa }+\stackrel{~}{D}_{ni}^{\lambda \mu }\pi _{in^{}}^\kappa )\left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right)+\frac{1}{2}\underset{i}{\overset{B}{}}\underset{j}{\overset{B}{}}\pi _{ni}^\lambda \pi _{ij}^\mu \pi _{jn^{}}^\kappa \left(\frac{1}{E_{ni}E_{nj}}+\frac{1}{E_{n^{}i}E_{n^{}j}}\right)\hfill \\ \hfill \frac{1}{2}\underset{i}{\overset{B}{}}\underset{n^{\prime \prime }}{\overset{A}{}}\left(\frac{\pi _{nn^{\prime \prime }}^\lambda \pi _{n^{\prime \prime }i}^\mu \pi _{in^{}}^\kappa }{E_{ni}E_{n^{\prime \prime }i}}+\frac{\pi _{ni}^\lambda \pi _{in^{\prime \prime }}^\mu \pi _{n^{\prime \prime }n^{}}^\kappa }{E_{n^{}i}E_{n^{\prime \prime }i}}\right),\end{array}$$
(115)
$$\begin{array}{c}Q_{nn^{}}^{\lambda \mu \kappa \nu }=\stackrel{~}{Q}_{nn^{}}^{\lambda \mu \kappa \nu }+\frac{1}{2}\underset{i}{\overset{B}{}}(\stackrel{~}{D}_{ni}^{\lambda \mu }\stackrel{~}{D}_{in^{}}^{\kappa \nu }+\pi _{ni}^\lambda \stackrel{~}{C}_{in^{}}^{\mu \kappa \nu }+\stackrel{~}{C}_{ni}^{\lambda \mu \kappa }\pi _{in^{}}^\nu )\left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right)\hfill \\ \hfill +\frac{1}{2}\underset{i}{\overset{B}{}}\underset{j}{\overset{B}{}}(\stackrel{~}{D}_{ni}^{\lambda \mu }\pi _{ij}^\kappa \pi _{jn^{}}^\nu +\pi _{ni}^\lambda \stackrel{~}{D}_{ij}^{\mu \kappa }\pi _{jn^{}}^\nu +\pi _{ni}^\lambda \pi _{ij}^\mu \stackrel{~}{D}_{jn^{}}^{\kappa \nu })\left(\frac{1}{E_{ni}E_{nj}}+\frac{1}{E_{n^{}i}E_{n^{}j}}\right)\\ \hfill \frac{1}{2}\underset{i}{\overset{B}{}}\underset{n^{\prime \prime }}{\overset{A}{}}\left(\frac{\stackrel{~}{D}_{nn^{\prime \prime }}^{\lambda \mu }\pi _{n^{\prime \prime }i}^\kappa \pi _{in^{}}^\nu +\pi _{nn^{\prime \prime }}^\lambda \stackrel{~}{D}_{n^{\prime \prime }i}^{\mu \kappa }\pi _{in^{}}^\nu +\pi _{nn^{\prime \prime }}^\lambda \pi _{n^{\prime \prime }i}^\mu \stackrel{~}{D}_{in^{}}^{\kappa \nu }}{E_{ni}E_{n^{\prime \prime }i}}+\frac{\stackrel{~}{D}_{ni}^{\lambda \mu }\pi _{in^{\prime \prime }}^\kappa \pi _{n^{\prime \prime }n^{}}^\nu +\pi _{ni}^\lambda \stackrel{~}{D}_{in^{\prime \prime }}^{\mu \kappa }\pi _{n^{\prime \prime }n^{}}^\nu +\pi _{ni}^\lambda \pi _{in^{\prime \prime }}^\mu \stackrel{~}{D}_{n^{\prime \prime }n^{}}^{\kappa \nu }}{E_{n^{}i}E_{n^{\prime \prime }i}}\right)\\ \hfill +O(\pi ^4).\end{array}$$
(116)
In the last expression, the symbol $`O(\pi ^4)`$ denotes terms of the fourth order in the kinetic momentum $`\pi `$. These cannot be obtained from Eq. (113), but they can be derived easily from Eq. (B15e) on p. 205 of Ref. Winkler, 2003.

The term $`W_{nn^{}}^{(2)}(𝐪)`$ is a renormalized second-order contribution to the band offsets, given by

$$W_{nn^{}}^{(2)}(𝐱)=\stackrel{~}{W}_{nn^{}}^{(2)}(𝐱)+\frac{1}{2}\underset{i}{\overset{B}{}}\left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right)[W_{ni}^{(1)}(𝐱)+\mathrm{\Lambda }_{ni}\phi ^{(1)}(𝐱)][W_{in^{}}^{(1)}(𝐱)+\mathrm{\Lambda }_{in^{}}\phi ^{(1)}(𝐱)],$$
(117)
where $`W_{nn^{}}^{(1)}(𝐪)`$ and $`\stackrel{~}{W}_{nn^{}}^{(2)}(𝐪)`$ were defined in Eqs. (64) and (68). The only contribution from $`\phi ^{(1)}`$ to be included here is the linear interface dipole term (67) for heterovalent systems; but note that this contribution is zero in the $`GW`$ approximation Hedin (1965); Hedin and Lundqvist (1969) and in density-functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965) (since $`\mathrm{\Lambda }_{nn^{}}=\delta _{nn^{}}`$ in these cases).

The remaining terms in (73) are renormalized versions of the functions $`\stackrel{~}{J}_{nn^{}}^\lambda (𝐪)`$, $`\stackrel{~}{M}_{nn^{}}^{\lambda \mu }(𝐪)`$, and $`\stackrel{~}{R}_{nn^{}}^{\lambda \mu }(𝐪)`$ defined in Eqs. (64) and (65):

$$J_{nn^{}}^\lambda (𝐪)=\stackrel{~}{J}_{nn^{}}^\lambda (𝐪)+\frac{1}{2}\underset{i}{\overset{B}{}}\pi _{ni}^\lambda W_{in^{}}^{(1)}(𝐪)\left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right),$$
(118)
$$\begin{array}{c}M_{nn^{}}^{\lambda \mu }(𝐪)=\stackrel{~}{M}_{nn^{}}^{\lambda \mu }(𝐪)+\frac{1}{2}\underset{i}{\overset{B}{}}\stackrel{~}{D}_{ni}^{\lambda \mu }W_{in^{}}^{(1)}(𝐪)\left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right)+\frac{1}{2}\underset{i}{\overset{B}{}}\pi _{ni}^\lambda \stackrel{~}{J}_{in^{}}^\mu (𝐪)\left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right)\hfill \\ \hfill +\frac{1}{2}\underset{i}{\overset{B}{}}\underset{j}{\overset{B}{}}\pi _{ni}^\lambda \pi _{ij}^\mu W_{jn^{}}^{(1)}(𝐪)\left(\frac{1}{E_{ni}E_{nj}}+\frac{1}{E_{n^{}i}E_{n^{}j}}\right)\\ \hfill \frac{1}{2}\underset{i}{\overset{B}{}}\underset{n^{\prime \prime }}{\overset{A}{}}\left(\frac{\pi _{nn^{\prime \prime }}^\lambda \pi _{n^{\prime \prime }i}^\mu W_{in^{}}^{(1)}(𝐪)}{E_{ni}E_{n^{\prime \prime }i}}+\frac{\pi _{ni}^\lambda \pi _{in^{\prime \prime }}^\mu W_{n^{\prime \prime }n^{}}^{(1)}(𝐪)}{E_{n^{}i}E_{n^{\prime \prime }i}}\right),\end{array}$$
(119)
$$\begin{array}{c}R_{nn^{}}^{\lambda \mu }(𝐪)=\stackrel{~}{R}_{nn^{}}^{\lambda \mu }(𝐪)+\frac{1}{2}\underset{i}{\overset{B}{}}\left(\pi _{ni}^\lambda [\stackrel{~}{J}_{n^{}i}^\mu (𝐪)]^{}+\stackrel{~}{J}_{ni}^\lambda (𝐪)\pi _{in^{}}^\mu \right)\left(\frac{1}{E_{ni}}+\frac{1}{E_{n^{}i}}\right)\hfill \\ \hfill +\frac{1}{2}\underset{i}{\overset{B}{}}\underset{j}{\overset{B}{}}\pi _{ni}^\lambda W_{ij}^{(1)}(𝐪)\pi _{jn^{}}^\mu \left(\frac{1}{E_{ni}E_{nj}}+\frac{1}{E_{n^{}i}E_{n^{}j}}\right)\\ \hfill \frac{1}{2}\underset{i}{\overset{B}{}}\underset{n^{\prime \prime }}{\overset{A}{}}\left(\frac{\pi _{nn^{\prime \prime }}^\lambda W_{n^{\prime \prime }i}^{(1)}(𝐪)\pi _{in^{}}^\mu }{E_{ni}E_{n^{\prime \prime }i}}+\frac{\pi _{ni}^\lambda W_{in^{\prime \prime }}^{(1)}(𝐪)\pi _{n^{\prime \prime }n^{}}^\mu }{E_{n^{}i}E_{n^{\prime \prime }i}}\right).\end{array}$$
(120)
Note that $`M_{nn^{}}^{\lambda \mu }(𝐪)M_{nn^{}}^{\mu \lambda }(𝐪)`$, in contrast to Eq. (66).