An envelope function formalism for lattice-matched heterostructures
Maarten L. Van de Put William G. Vandenberghe Wim Magnus Bart Sorée
cond-mat.mes-hall




.eps


    


ua,imec]Maarten Van de Put
maarten.vandeput@imec.be
utd]William Vandenberghe
ua,imec]Wim Magnus
ua,imec]Bart Sorée

[ua]Department of Physics, Universiteit Antwerpen, Antwerpen, Belgium
[imec]imec, Kapeldreef 75, 3001 Heverlee, Belgium
[utd]Department of Materials Science and Engineering, The University of Texas at Dallas, Dallas, USA



The envelope function method traditionally employs a single basis set which, in practice, relates to a single material because the · matrix elements are generally only known in a particular basis. 
In this work, we defined a basis function transformation to alleviate this restriction. 
The transformation is completely described by the known inter-band momentum matrix elements. 
The resulting envelope function equation can solve the electronic structure in lattice matched heterostructures without resorting to boundary conditions at the interface between materials, while all unit-cell averaged observables can be calculated as with the standard envelope function formalism. 
In the case of two coupled bands, this heterostructure formalism is equivalent to the standard formalism while taking position dependent matrix elements.



    ·, heterostructures, envelope functions
    73.22.-f






§ INTRODUCTION


To calculate the band structure of semiconductors, a number of methods
are available: the k· p method<cit.>,
the tight-binding method<cit.> and the pseudopotential
method<cit.>. The latter two methods are
atomistic while the k· p method is not and will
be computationally more efficient when atomic resolution is not required.
The k· p method has been successfully used to explain
many physical phenomena such as band-to-band tunneling (BTBT)<cit.>
or topological insulators<cit.>. The k· p
method has also been used successfully to study tunnel field-effect
transistors (TFETs)<cit.>. But with a recent interest
in TFETs composed of different materials<cit.>,
the need for a k· p-based framework that can deal
with heterostructures emerges.

Historically, the · formalism for bulk has been extended in various ways to an envelope function formalism which can provide the electronic structure for arbitrary potentials.
The rigorous derivation of an exact envelope function formalism by Burt<cit.> provides a solid physical and mathematical basis for most extensions of the · formalism.
The exact envelope function equations described in Burt's work involve a non-local interaction of the external potential, but under the right approximations (a slowly varying potential) they reproduce local equations. 
In literature the local equations are used almost exclusively because they are numerically manageable due to known, empirically determined, matrix elements while Burt's exact equations require knowledge of the complete basis set.

Several authors have proposed strategies to adopt the local envelope function formalism to heterostructures, we recognize three common approaches.
The first approach correctly accounts for the different material basis sets by relying on a complete knowledge of the basis functions<cit.>. 
This method is a nice theoretical exercise, but highly impractical, as most benefits of the envelope function are lost when full knowledge of the basis functions is required.
A second strategy (implicitly) assumes the basis functions for different materials used for the envelope function expansion to be identical for different materials, which enables bulk-like envelope function expansion. However, in practice the matrix elements are often taken to be the local material specific bulk matrix elements in their natural basis<cit.>, which is inconsistent with them being identical.
Foreman details a third strategy, he accounts for the difference in basis functions by introducing a interface term to describe additional couplings between the two materials, while using material-dependent basis functions throughout.<cit.> This last strategy has recently been employed in the context of studying the homogenization limit for heterostructure · multi-band models and the derivation of an optimal effective mass model for heterostructures<cit.>.

In this paper we introduce a new set of envelope function equations which is capable of describing the electronic structure in lattice-matched heterostructures, taking full account of the difference in basis sets. 
The materials can have very different bulk electronic structures, which are properly handled, equivalent to the Foreman technique.
However, we do so by expanding on a single basis set, which leads to a continuous set of equations free of interface terms, which seems appropriate given a single Hilbert space for the solution to the entire heterostructure.
Neither new parameters, nor additional assumptions regarding the material properties are in order,
only the standard · matrix elements of every material involved need to be known.
A benchmark using this new formalism to calculate band-to-band transitions (BTBT) in heterostructures has been presented and published<cit.>, and more recently, a two dimensional BTBT simulator has been implemented<cit.>.
In the present paper, we elaborate on the underlying theory and its implications.



§ HETEROSTRUCTURE MODEL: THE SCHRÖDINGER EQUATION


The goal of this paper is to develop an envelope function formalism to solve the one-electron Schrödinger equation in lattice-matched semiconductor heterostructures under the Hartree approximation,

    -ħ^2/2m_0^2 ψ() + V_ c() ψ() + V_ e() ψ() = E ψ(),

where the crystal potential V_ c() is taken piecewise throughout the structure and the extrinsic potential V_ e() contains all contributions to the potential energy not related to the bulk material.

We distinguish between the heterostructure materials with an index λ, each material having a region Ω_λ in the heterostructure and a boundary ∂Ω_λ as depicted in Fig. <ref>.


Using this notation we can write the full form of the crystal potential,

    V_ c() = ∑_λθ(∈Ω_λ) V^λ_ c(),

where the sum runs over all materials λ. θ(∈Ω_λ) is the logical step function, or indicator function, which equals 1 if  is in Ω_λ, 0 outside this region and it selects the correct crystal potential for a given region λ. The crystal potentials V_ c^λ() are lattice periodic and defined in all space.



§ AN ARBITRARY BASIS

We expand the wavefunction in an arbitrary set of complete, periodic basis functions u_n(), yielding

    ψ() = ∑_n f_n() u_n(),

where the expansion coefficients f_n() are called envelope functions. Due to the periodicity of the basis functions, this expansion is unique if the plane wave expansion of the envelope functions only has non-zero components within the first Brillouin zone.

Substitution into the Schrödinger equation Eq. (<ref>) leads to

    ∑_n [ 
        -ħ^2/2m_0^2 f_n() u_n()
        - ħ/m_0 f_n() · (-ħ) u_n() 
        + f_n() -ħ^2/2m_0^2 u_n() 
    
        + V_ c() f_n() u_n() 
        + V_ e() f_n() u_n()
      ]
      = ∑_n E  f_n() u_n().

In the next sections we introduce momentum and energy matrix elements in order to recast Eq. (<ref>) into a set of coupled equations in which the differential operators merely act on the envelope functions.



 §.§ Momentum term

Starting with the second term of Eq. (<ref>), we write -ħ u_n() as an integral over a unit cell,

    -ħ u_n() = ∫_ uc^3r' δ(-') (-ħ') u_n('),

and exploit the completeness relation ∑_m u^*_m(') u_m() = δ(-'), yielding

    -ħ u_n() = ∑_m ∫_ uc^3r' 
          u_m(') (-ħ') u_n(') u_m().

where the periodicity of the basis functions restricts the integration to a single unit cell (uc).
Substitution into the second term of Eq. (<ref>) and swapping n and m leads to

    ∑_nm -ħ/m_0 f_m() ·_nm u_n().

The momentum matrix elements are empirically determined for basis functions which are solutions for a bulk material corresponding to a point of high symmetry in the first Brillouin zone,

    _nm = -ħ∫_ uc d^3r  u^*_n()  u_m().




 §.§ Bulk Hamiltonian term


Combining the third an fourth term of Eq.(<ref>), we recover the bulk Hamiltonian at ǩ=(000) of layer λ,

    ∑_n [ -ħ^2/2m_0^2 + V_ c() ] u_n() f_n()

Expanding the crystal potential as in Eq. (<ref>) and using the fact that ∑_λθ(∈Ω_λ) = 1 we obtain

    ∑_n ∑_λθ(∈Ω_λ) 
        [ -ħ^2/2m_0^2 + V^λ_ c() ] u_n() f_n()

Inserting a completeness relation for the two terms and swapping indices n and m, we obtain an expression in the basis functions directly,

    ∑_nm∑_λθ(∈Ω_λ) H^λ_nm f_m()
      u_n(),

where the matrix elements of the bulk Hamiltonian for layer λ are given by

    H^λ_nm = ∫_ uc^3r  u^*_n() [ -ħ^2/2m_0^2 + V^λ_ c() ] u_m().




 §.§ Envelope function equation

Using the completeness of the u_n(), we arrive at a set of envelope function equations

    - ħ^2/2m_0^2 f_n()
      - ħ/m_0∑_m _nm· f_m() 
      + ∑_m H_nm() f_m() 
      + V_ e() f_n() 
      = E  f_n(),

with momentum matrix elements _nm as defined in Eq. (<ref>), and the position dependent bulk-Hamiltonian matrix elements

    H_nm() = ∑_λθ(∈Ω_λ) H^λ_nm,

which use the material dependent bulk Hamiltonian matrix elements H^λ_nm as in Eq. (<ref>).

Because the envelope functions do not have to be restricted to the first Brillouin zone, the expansion of the wavefunctions as specified in Eq. (<ref>) is not unique, but rather convenient. A derivation of envelope functions which are unique, but non-local by construction, is detailed in the appendix.



§  SELECTING A BASIS SET: THE HETEROSTRUCTURE PROBLEM

As mentioned before, we need to select a specific set of basis functions such that the matrix elements can be empirically determined. 

Choosing the u_n^ρ() basis to be the zone-centre solutions of the bulk material denoted by index ρ we may write the envelope function expansion of the wavefunction as

    ψ() = ∑_n f^ρ_n() u^ρ_n(),

whereas the basis dependent envelope functions f^ρ_n() satisfy the envelope function equations given by

    - ħ^2/2m_0^2 f^ρ_n()
      - ħ/m_0∑_m ^ρ_nm· f^ρ_m() 
      + ∑_m H^ρ_nm() f^ρ_m() 
      + V_ e() f^ρ_n() 
      = E  f^ρ_n().

Here, the inter-band momentum matrix elements ^ρ_nm only depend on the chosen basis and are known for many basis sets.<cit.> Being independent on the local material, they are position independent.

The bulk Hamiltonian matrix elements are dependent on the basis through ρ, and the local crystal potential through λ,

    H_nm^ρ,λ = ∫^3r  u^ρ*_n() 
        [ -ħ^2/2m_0^2 + V^λ_ c() ] u^ρ_m().

They are well known when ρ=λ, namely

    H^λ,λ_nm = δ_n,m E^λ_n,

where E^λ_n are the band zone-centre energies associated with the Bloch basis {u_n^λ()}.

However, when considering H_nm^ρ,λ with ρ≠λ we cannot rely on empirically obtained data and one would need to evaluate the integral in Eq. (<ref>), as this would require full knowledge of the zone-centre basis set {u_n^ρ()} and the mean-field crystal potentials V^λ_ c() for all λ. Moreover the integral would have to be computed numerically in general.



§ BASIS FUNCTION TRANSFORMATION


To circumvent the problem of the unknown matrix elements we perform a unitary transformation S^λ→ρ which transforms the u_m^λ() basis into u_n^ρ(), i.e.

    u^ρ_n() = ∑_m S_nm^λ→ρ  u^λ_m(),

where

    S_nm^λ→ρ = ∫ d^3r   u^ρ*_n() u^λ_m().


We could, in principle at least, determine the transformation matrix elements from the set of basis functions {u^ρ_nm()} by numerical integration of Eq. (<ref>). However, a central point to the practical use of · theory is that we don't need to know the basis functions. Finding the basis function requires the solution of the full bulk Schrödinger equation at the =(000) point, which in turn requires knowledge of the crystal potential V_ c^λ() ∀λ. With ·, we avoid this problem by introducing empirically determined inter-band momentum matrix elements p_nm and zone-centre energies E_n. In this section we will detail a method in which the elements of the transformation matrix S_nm^λ→ρ are be obtained solely from these empirical parameters.



 §.§ The transformation matrix elements


We set out by introducing a basis set of plane wave with lattice periodicity, namely 

    u_() = e^i·,

where  is a reciprocal lattice vector.

As the plane wave basis functions are eigenfunctions of the momentum operator the inter-band momentum matrix elements turn out to be diagonal,

    _' 
      = -ħ∫ d^3r  u^*_()  u_() 
      = δ_' ħ.


Eigenvalue decomposition of the inter-band momentum matrix P̌^λ of any complete set of basis functions u_n^λ() with the same lattice periodicity will produce these exact same eigenvalues:

    P̌^λ = R^λ  diag(ħ) R^λ†,

where we assume the eigenvectors in the columns of R^λ are properly normalized so that R^λ is a unitary matrix. 
The momentum matrix elements read

    _nm^λ 
      = ∑_ R_n^λ ħ  R_ m^λ†
      = ∑_ R_n^λ ħ  R_m^λ*

where the matrix elements R_n^λ are the eigenvectors of the inter-band momentum matrix P̌^λ. 

By defining a transformation Q^λ from the plane wave basis u_() to the basis functions u^λ_n() of material λ,

    u^λ_n() = ∑_ Q^λ_nG u_(),

and using this transformation on the inter-band momentum matrix, we obtain a similar equation as in Eq. (<ref>),

    ^λ_nm   = -ħ∫ d^3r  u^λ*_n()  u^λ_m() 
       = -ħ∑_'∫ d^3r  
        Q^λ*_n u^*_()  Q^λ_m u_'() 
       = ∑_ Q^λ*_n ħ Q^λ_m,

or in matrix-representation

    P̌^λ = Q^λ* ħ Q^λ⊺,

from which we conclude R^λ=Q^λ*. The matrix elements Q^λ_n of the transformation from plane waves to the basis set of material λ are completely determined by the eigendecomposition of the inter-band momentum matrix of the corresponding material. The inter-band momentum matrix elements can be obtained from · models with a high number of bands, readily available in literature for common materials.<cit.>
  
For the transformation between two material basis sets u^λ_n() and u^ρ_n(), we derive

    S_nm^λ→ρ   = ∫^3r  u_n^ρ*() u_m^λ() 
      = ∫^3r ∑_∑_' 
        Q_n^ρ* e^-· Q_m'^λ e^·
       = ∑_∑_' Q_n^ρ*  Q_m'^λδ_,' 
      = ∑_ Q_n^ρ*  Q_m^λ

or in matrix format:

    S^λ→ρ = Q^ρ†  Q^λ.


In order to obtain the transformation matrix S^λ→ρ we do not need to evaluate the integral in Eq. (<ref>). Instead, we need only to diagonalize the inter-band momentum matrix P̌^λ so as to obtain Q^λ, while the inter-band momentum matrix is a well known object. 



 §.§ Revisiting the bulk Hamiltonian matrix elements


With the transformation matrix elements known, we can now easily express the transformation matrix elements H^λ,ρ_nm of Eq. <ref> in terms of the available material specific matrix elements by transforming from the reference basis set to the zone-centre solutions u^λ_n(),

    H^λ,ρ_nm   = ∫ d^3r  
          u^ρ*_n() [-ħ^2/2m^2 + V_ c^λ() ] u^ρ_m()
       = ∑_ij∫ d^3r  
          u^λ*_i() [S^λ→ρ_ni]^*
          Ĥ^λ S^λ→ρ_mj u^λ_j()
       = ∑_ij[S^λ→ρ]^*_ni H^λ,λ_ij S^λ→ρ_mj

Because the bulk Hamiltonian of material λ is now expressed in it's own natural basis set the matrix elements reduce to the band energies at the zone-centre H^λ,λ_ij = δ_ij E^λ_i. The bulk-Hamiltonian matrix elements of material λ in the u^ρ_n() basis are thus found as

    H^λ,ρ_nm = ∑_i [S^λ→ρ_ni]^* E^λ_i S^λ→ρ_mi

As expected, the columns of the complex conjugated transformation matrix are the eigenvectors of the bulk-Hamiltonian of a material λ expressed in the reference basis u_n^ρ(). This is equivalent to the momentum matrix decomposition that led to the transformation matrix elements.



 §.§ Transforming the envelope functions


The envelope functions are expansion coefficients with respect to the chosen reference basis set. In most applications, it is however convenient to have envelope functions defined against the material's own zone-centre bulk-Hamiltonian solutions. We achieve this by transforming the calculated envelope functions in each material to it's proper material specific basis set.

In the region of material λ we transform to the basis functions u^λ_n(),

    ψ() = ∑_m f^ρ_m() u^ρ_m() 
                = ∑_mn f^ρ_m() S^λ→ρ_mn u^λ_n()
                = ∑_n f^λ_n() u^λ_n(),

which amounts to

    f^λ_n() = ∑_m f^ρ_m() S^λ→ρ_mn.

These envelope functions can be interpreted just like those produced by the homostructure envelope function formalism.



§ BASIS CHOICE AND RESTRICTION


While the basis transformation enables envelope function expansion with respect to a particular reference basis set, one is free to select that reference basis set.
In this section, we investigate the effects of this choice when the set of basis functions is truncated such that the number of basis functions is finite. 
The truncated basis set is no longer complete and does not span the whole solution space of the original Schrödinger equation.
When changing basis sets, we have two options for the momentum matrix elements: either the momentum matrix elements are transformed like the Hamiltonian matrix elements, i.e. P^λ→[S^λ→ρ]^† P^λS^λ→ρ, or the momentum matrix elements P^ρ from the new basis are used.
We show that the former option leads to a set of equations that is invariant under basis choice.



 §.§ A two band model




In this section, we study the simplest case by restricting the basis set to just two functions u^ρ_1() and u^ρ_2(), exhibiting a non-zero inter-band momentum matrix element between them only in the x direction,

    ^ρ_12 = -ħ∫^3r  u_1^ρ*() ∇ u^ρ_2() = (p^ρ,0,0).

One can think of the conduction band state with s-like symmetry and a valence band state with p_x-like symmetry in the Γ-point of semiconductors. In this simple model a major problem arises due to the freedom of choice for the reference basis set. This problem will be resolved in the next section

The coupling being one-dimensional, Eq. (<ref>) now reduces to

    - ħ^2/2m_0x
      - ħ/m_0 P^ρx
      + H^ρ(x) 
      + V_ e(x) 
      = E ,

where the envelope function are f_1^ρ(x) and f_2^ρ(x), with a simple inter-band momentum matrix,

    P^ρ = [[    0  p^ρ; p^ρ*    0 ]].


In this oversimplified case it is not even necessary to calculate the transformation matrix elements as the symmetry considerations lead exclusively to the identity transformation,

    S^ρ→λ = [[ 1 0; 0 1 ]],

for every lattice matched material combination λ, ρ. 

In the two-band case, when using our basis transformation model, there is no distinction between the basis functions of the different materials and this results in a very simple set of envelope function equations. Only the basis function energies E_1 and E_2 are material dependent,

    {- ħ^2/2m_0x
      - ħ/m_0[[       0     p^ρ; (p^ρ)^*       0 ]] x 
      + ∑_λθ(x∈Ω_λ) 
          [[ E_ 1^λ      0;      0 E_ 2^λ ]]
      + V_ e(x) }
      = E .

The inability of the transformation to distinguish the basis functions of two different materials and the envelope function equations defined earlier has introduced an ambiguity. Changing the basis changes the inter-band momentum matrix elements but not the band edge energies and results in different equations with different solutions.

To make the ambiguity explicit, we consider a bulk material λ, without any external potential. The envelope function equations in the basis set u^λ_n(x) yield the expected · equations,

    {- ħ^2/2m_0x 
      - ħ/m_0[[       0     p^λ; (p^λ)^*       0 ]] x
      + [[ E_ c^λ      0;      0 E_ v^λ ]]
      + V_ e(x) }
      = E

but transforming the into a different basis set u_n^ρ(x), Eq. (<ref>) becomes

    {- ħ^2/2m_0x
      - ħ/m_0[[       0     p^ρ; (p^ρ)^*       0 ]] x
      + [[ E_ c^λ      0;      0 E_ v^λ ]]
      + V_ e(x) }
      = E

and in general p^λ≠ p^ρ, so our model seems to provide the wrong · equation for bulk material λ when expressed in any basis set u_n^λ(x) other than its own. 

Hence, to remedy this problem, we should use a material dependent inter-band momentum matrix element, giving rise to the following two band envelope function equations:

    {- ħ^2/2m_0x
      - iħ/m_0∑_λθ(x∈Ω_λ)
        [[       0     p^λ; (p^λ)^*       0 ]] -ħ/m_0x
      + ∑_λθ(x∈Ω_λ)
        [[ E_ c^λ      0;      0 E_ v^λ ]] 
         
      + V_ e(x) }
    
      = E .

One can easily check that these equations do produce the right · equations for bulk material, for any basis.

In the two band case, our heterostructure formalism with basis function transformations is equivalent to taking position dependent matrix elements in the classical envelope function method without any transformations.



 §.§ Robust envelope function equations

The above considerations clearly indicate that using a finite number of basis functions calls for transformed momentum matrix elements in the envelope function equations. In the infinite complete basis set we have

    ^ρ_nm = ∑_ij[S^λ→ρ_ni]^* ^λ_ij S^λ→ρ_mj.

with which we obtain a modified set of envelope function equations,

    - ħ^2/2m_0^2 f^ρ_n()
      - ħ/m_0∑_m ^ρ_nm() · f^ρ_m() 
      + ∑_m H^ρ_nm() f^ρ_m() 
      + V_ e() f^ρ_n()
      = E  f^ρ_n(),

with position dependent momentum matrix elements,

    ^ρ_nm()    = ∑_λθ(∈Ω_λ)  ^λ,ρ_nm,
    ^λ,ρ_nm   = ∑_ij[S^λ→ρ_in]^* p^λ_ij S^λ→ρ_jm.

While, in the case of a complete basis set, the exact relation ^λ,ρ_nm = ^ρ_nm still holds, using this new, robust form of envelope equation the choice of reference basis set has no influence on the final results, even when these basis sets are as heavily restricted as in the two band model.
This can be verified for any arbitrary restriction by performing a transformation of the envelope functions f_n^ρ() to the f_m^ρ() basis set which results in pairs of transformation matrices [S^λ→ρ]^† S^λ→ρ that cancel because of unitarity, even for an incomplete basis set.

Furthermore, because of the symmetry constraints imposed on the basis set, this form results in the same equations used in literature whenever the set is restricted to bands with all different symmetry as is the case in the two band model, but also in the four (eight) band model of a semiconductor with conduction band states with s-like symmetry and valence states  with p_x-, p_y- and p_z-like symmetry. However, for more accurate models such as those including 15 or 30 bands<cit.>, this form will provide more accurate results, certainly when inter-band processes are important.



§ UNIT-CELL OBSERVABLES


Being represented by operators acting on the entire wavefunction, observables are independent of the chosen basis.
It is however desirable to define operators acting only on the envelope function. These can only give approximate results for the observables because they are defined by averaging over one unit-cell to remove the basis function dependence.
In this section, we show expressions for the unit-cell average density and probability current.
We show that these observables can be determined from the envelope functions in any basis set without explicit transformation to the original material basis.

We define the average unit-cell density of a pure state ψ() as

    <ρ>_ uc()=∑ _nf_n^ρ*()f_n^ρ()

and since any transformation to another basis set f_n^ρ() is unitary, the average unit-cell density is basis independent.

Next, we turn to the probability current. The probability current is given by

    J̌_ prob() = -ħ/2m_0[ ψ^*() ψ() - ψ() ψ^*() ].

We use the envelope function expansion and take an average over one unit cell, removing the basis function dependence and introducing inter-band momentum matrix elements in the equation,

    <J̌_ prob>_ uc() =
        - ħ/2m_0∑_n 
          [ f^ρ*_n()  f^ρ_n() - f^ρ_n()  f^ρ*_n() ] 
        + 1/m_0∑_n,m^ρ_nm(r) f^ρ*_n() f^ρ_m().

Here we have again used the position dependent matrix element to make Eq. (<ref>) robust to basis changes when considering an incomplete basis set. Upon transformation of the envelope functions to another basis we find that it is basis independent due to unitarity of the transformation.



§ CONCLUSIONS

We have proposed a transformation between zone-centre solutions of the bulk Hamiltonian of different materials. The transformation is fully described by the eigendecomposition of the inter-band momentum matrices of both materials. The inter-band momentum matrices are empirically determined by · fitting to bulk material and available in literature for many materials.

This transformation enables the expansion of the wave function on the zone-centre solutions of any particular material. 
Without this transformation we are confronted with unknown matrix elements that cannot be empirically determined. 
In our approach, we determined the unknown matrix elements by transforming them from a natural basis set where they are known to a singular reference basis set which is used as a basis for the envelope expansion. No new parameters were needed in this process.

Our method performs a similar function to previously published work, and differs from it mostly in a practical sense. For example the original formalism of Burt<cit.> allows for heterostructures and takes account of non-local effects due to abrupt changes in potential. However, in application of this method it is assumed that the basis functions are known. As stated earlier this is not the case for empirically fitted · models. Our method is also closely related to the basis-dependent method by Foreman<cit.> as used in <cit.>, where the basis functions change with the material while interface terms mediate the coupling (transformation) from one material to the next. Our method effectively uses a single basis set throughout the whole structure and instead transforms the · parameters in each region to this shared basis. Because we expand on a single basis set no explicit interface terms appear.

We should note that the non-local caused by abrupt changes of the potential present in the formalisms by Burt and Foreman are also applicable to our method and can be found in the Appendix. This also includes non-local effects caused by abrupt changes in momentum and energetic coupling between materials.

We have taken care to make the heterostructure envelope functions robust under change of basis when this basis is no longer complete. This method is thus applicable even for only a few bands with the same accuracy as the homostructure based approach. For the basis set of two coupled bands, and other basis sets with heavy symmetry restrictions, we obtained a simpler set of equations, where in each region the homostructure equation is retrieved and no transformation is present.

To conclude, we note that the computational burden related to this method is no greater than for the conventional homostructure methods. The transformation matrix calculation involves a computationally inexpensive simultaneous eigenvalue decomposition of the 3D inter-band momentum matrix. The transformation matrix can even be tabulated in advance of any calculations.



§ APPENDIX: UNIQUE ENVELOPE FUNCTIONS



The heterostructure envelope function equations discussed in this paper do not correspond to a unique expansion of the wavefunction. As we show in this section, a unique expansion results in non-local integro-differential envelope function equations, which are not as convenient. The convenient, local, if not unique, equations discussed throughout the paper are based on the slowly-varying field approximation commonly used in practical envelope function calculations. This means that the external potential and the matrix elements must not change abruptly on the scale of a unit cell. For the external potential, this should be evaluated on a case-by-case basis. As for the matrix elements, we can appeal to physical demands that the basis set and crystal potential do not contain abrupt discontinuities due to the non-local nature of the electron interaction with the lattice.

It is however possible to derive envelope function equations for a unique expansion of the basis set. We enforce this uniqueness by restricting the non-zero plane wave components of the envelope functions to the first Brillouin Zone (1BZ). 
 We obtain envelope function equations that return exact and unique envelope functions at the cost of locality of the envelope function equations.



 §.§ External potential

For the external potential, we adopt Burt's methodology to restrict the plane wave components to the 1BZ, yielding<cit.>

    ∑_n V_ e() f_n() u_n() 
    = ∑_nm∫^3r' V_nm(,') f_m(') u_n(),

where the non-local kernel V_nm(,') is

    V_nm(,') = ∑_'∑_'[ ũ_n+_1]^* Ṽ_-'() ũ_m' 
    e^(_1·-'·').

Here, ũ_n and Ṽ_() are the Fourier transform of u_n() and V_ e() respectively. _1 is a wavevector inside the 1BZ and _1 is a reciprocal lattice vector defined by +'=_1+_1.

Invoking a basis set transformation, we have now also determined the Fourier components by the eigendecomposition of the momentum matrix as ũ^ρ_n = Q_n. The exact knowledge of the basis functions is thus no longer needed to determine the non-local interaction.



 §.§ Bulk Hamiltonian matrix elements

We apply a similar procedure to the bulk Hamiltonian matrix elements, starting with the expansion of Eq. (<ref>) into plane waves,

    ∑_nm H_nm() f_m() u_n()
    = ∑_nm∑_'∑_'H̃_nm,(k) f̃_m(') ũ_n' e^(++'+')·.

Substituting →-', expressing +' as a 1BZ restricted wavevector _1, and introducing a reciprocal lattice vector _1 with +=_1+_1,

    ∑_nm∑_'∑_GG'H̃_nm,-'(k) f̃_m(') ũ_n' e^_1· e^(+_1)·

we finally rewrite e^(+_1)· using the basis functions e^(+_1)· = ∑_j ũ_j,+_1^* u_j(),

    ∑_nmj∑_'∑_'H̃_nm,-'() f̃_m(') ũ_n'ũ_j,+_1^* e^_1· u_j().

We obtained an expansion on the basis functions u_j() where the expansion coefficients are restricted to the 1BZ because _1 is always inside the 1BZ. Performing a inverse Fourier decomposition on the envelope functions f̃_m(') = ∫^3r' f_m() e^-'·' and exchanging indices n and j yields

    ∑_nm H_nm() f_m() u_n()
    = ∑_n [ ∑_m ∫^3r' H_nm(,') f_m(') ] u_n()

with the position dependent, non-local Hamiltonian matrix elements,

    H_nm(,') = ∑_j ∑_'∑_'ũ_n,+_1^* H̃_jm,-'() ũ_j' e^(_1·-'·')
 
Here again we can use the known transformation matrix elements Q_n=ũ_n.



 §.§ Envelope functions

With Eq. (<ref>) and Eq. (<ref>) we obtain a unique expansion in terms of u_n(). The envelope functions are then given by a set of non-local envelope functions,

    - ħ^2/2m_0^2 f_n()
      - ħ/m_0∑_m _nm· f_m() 
      + ∑_m ∫^3r'[ H_nm(,') + V_nm(,') ] f_m(')
      = E  f_n().

These equations are much harder to solve due to the non-local interactions and the extra coupling due to the external potential. However, using the eigendecomposition detailed in this paper, the needed elements are all known and the equations are at least solvable in principle. 




§ REFERENCES

elsarticle-num


