Numerical solution of multiband model for tunnelling in type-II heterostructures

A new and very general method was developed for calculating the charge and spin-resolved electron tunnelling in type-II heterojunctions. Starting from a multiband description of the bulk energy-band structure, a multiband Riccati equation was derived. The reflection and transmission coefficients were obtained for each channel by integrating the Riccati equation over the entire heterostructure. Numerical instability was reduced through this method, in which the multichannel log-derivative of the envelope function matrix, rather than the envelope function itself, was propagated. As an example, a six-band Hamiltonian was used to calculate the current-voltage characteristics of a 10-nm wide InAs/ GaSb/InAs single quantum well device which exhibited negative differential resistance at room temperature. The calculated current as a function of applied (bias) voltage was found to be in semiquantitative agreement with the experiment, a result which indicated that inelastic transport mechanisms do not contribute significantly to the valley currents measured in this particular device. : spin-dependent electron transport, semiconductor nanostructures


Introduction
Electron tunnelling phenomena associated with the brokengap 1 configuration in type-II semiconductor heterostructures have been studied both theoretically 2 and experimentally. 3In these structures, the minimum in the conduction band of the material on one side of the hetero-interface lies at a lower energy than the maximum in the valence band of the material on the other side.Charge transfer between two adjacent layers can lead to the formation of a 2-D electron gas on one side and a 2-D hole gas on the other, even at zero bias.Interband tunnelling may thus occur between the spatially-separated intrinsic carriers situated on opposite sides of the hetero-interface, i.e. between electrons in the conduction band on one side and holes in the higher-lying valence bands on the other side. 4Under suitable conditions, electron-hole bound states (excitons) may be formed.The InAs/GaSb system, for example, is currently the most promising candidate for observing Bose-Einstein condensation of excitons. 5ince the electron states of both conduction and valence bands of the constituent semiconductors participate in the tunnelling process on an equal footing, theories of electron tunnelling in type-II heterostructures must include the band-mixing effects. 6Direct numerical solution of the multiband k•p Schrödinger equation, which is frequently employed to model these systems, is not feasible due to the inherent numerical instability that occurs in type-II systems.The numerical instability in these systems arises because of the simultaneous presence of propagating (real wave vector) and evanescent (imaginary wave vector) states.To circumvent the numerical difficulties that arise in such calculations, several innovative techniques have been proposed within the context of multiband k•p theory. 7,8ultiband formulations of the transfer matrix method (TMM) become numerically unstable in the simultaneous presence of evanescent and propagating states.For this reason the TMM is not very useful for studying type-II systems.In cases where the TMM fails, it is possible to propagate the scattering matrix instead. 9This technique was first applied to resonant tunnelling in GaAs/A1 x Ga 1-x As multilayer systems, with the influence of higher bands included. 10Although the scattering matrix method is numerically stable it is computationally no more efficient than the TMM.A third method, the multiband quantum transmitting boundary method, 11 is slightly more efficient; although with modern processor speeds the time saved on a typical calculation may be minimal in practice.
In Section 2.1 of this article, a new and very general method is proposed for calculating the electronic structure and transport properties of type-II heterostructures.In common with the preceding methods, the theory is based on a realistic, multiband k•p envelope function description of the energy-band structure.Starting from the general form of the n×n matrix Schrödinger equation, a multiband k•p Riccati equation is derived for the logarithmic derivative of the envelope function matrix.In Section 2.2 a boundary condition is derived for the Riccati equation.This boundary condition is used to integrate the Riccati equation numerically over the entire heterostructure.In Section 2.3 it is shown how to obtain the reflection matrix R, which contains the tunnelling amplitudes for each of the n channels along its diagonal.In Section 2.4 it is shown how the solution to the Riccati equation may be incorporated into a self-consistent calculation of the confining potential profile.The current through a device may then be calculated according to the definition given in Section 2.5.In Section 3, the room temperature current was calculated as a function of applied voltage for a 10-nm wide InAs/GaSb/InAs quantum well device and the results were compared to an experiment. 3It was found that the calculated current-voltage characteristics were in semi-quantitative agreement with the measured characteristics.This agreement suggests that, contrary to previous claims, inelastic transport mechanisms did not play a significant role in the resonant interband coupling mechanism that was most likely responsible for the negative differential resistance.

Multiband k•p Riccati equation
Consider a heterostructure with its growth direction along the x-axis.In multiband k•p models the energy-band structure within each layer of the heterostructure is described by a set of coupled differential equations for the components of the slowly varying envelope function, Φ(x).It can be shown that the Schrödinger equation for the envelope function can be decomposed into the following very general matrix form: 12 In Equation 1 the eigenvector Φ(x) is a column vector of length equal to the number of bands (n) which are explicitly included in the k•p model being used.The n × n matrices H 0 , H 1 and H 2 are usually derived through second order quasi-degenerate perturbation theory (Löwdin partitioning). 13Examples of the decomposition in Equation 1 can be found in Liu et al., 7 Harrison 8 and Eppenga et al. 14 Department of Physics, University of South Africa, P.O.Box 392, UNISA 0003, South Africa.E-mail: bothaae@unisa.ac.za

A new and very general method was developed for calculating the charge and spin-resolved electron tunnelling in type-II heterojunctions. Starting from a multiband description of the bulk energy-band structure, a multiband
Riccati equation was derived.The reflection and transmission coefficients were obtained for each channel by integrating the Riccati equation over the entire heterostructure.Numerical instability was reduced through this method, in which the multichannel log-derivative of the envelope function matrix, rather than the envelope function itself, was propagated.As an example, a six-band Hamiltonian was used to calculate the current-voltage characteristics of a 10-nm wide InAs/ GaSb/InAs single quantum well device which exhibited negative differential resistance at room temperature.The calculated current as a function of applied (bias) voltage was found to be in semiquantitative agreement with the experiment, a result which indicated that inelastic transport mechanisms do not contribute significantly to the valley currents measured in this particular device.
: spin-dependent electron transport, semiconductor nanostructures Instead of working with an n-component column vector Φ(x), it is advantageous to construct an n × n envelope function matrix, denoted by F(x).This matrix contains in its columns the n, linearly independent solution vectors of the stationary problem posed by Equation 1.
It can be shown quite generally that H 2 is non-singular. 15quation 1 can therefore be multiplied throughout by  5, Y will be referred to as the log-derivative of F.

Boundary condition for the Riccati equation
Consider the multiband potential profile of the heterostructure shown schematically in Fig. 1.The device consists of a central active region sandwiched between two flat-band contact regions.Assume that within the active region (x L , x R ) the potential may vary arbitrarily, but that it remains constant (not necessarily zero) in the outer two contact regions, i.e. for x ≤ x L and x ≥ x R .By making the very reasonable assumption that all charge carriers are incident from the left contact region, it follows that there can be no reflected waves in the right contact region.For (x ≥ x R ), the matrix solution to Equation 1 can therefore be written in terms of transmitted plane waves as Here D + , T and C are n × n matrices.D + is a known matrix whose columns are made up of the n linearly-independent eigenvectors of Equation 1 within the right contact region.Although D + can be obtained analytically for relatively small values of n (for example, see Botha 4 ), a more general numerical method has recently been suggested by Harrison. 8The numerical method consists of re-writing the non-linear eigenproblem posed by Equation 1 as a linear eigenproblem for the wave vector components (k i ) x (i = 1, 2,..., n) (see Equation 10.35 in Harrison 8 ).In Equation 6, T is an (as yet) unknown transmission matrix which depends only on the electron energy E, the transverse component of the carrier wave vector k || and the applied voltage V a .The matrix C is constant with respect to x and may be chosen arbitrarily.It is most convenient to choose C such that F R (x 0 ) = 1, where 1 is the n × n identity matrix and x 0 > x R .By making this choice, and evaluating Equation 6 at x = x 0 , it follows that Because the wave vector components (k i ) x and matrices D ± only depend on x indirectly through the potentials V i (x), D ± can be treated as constants within the two contact regions.Equation 6can therefore be differentiated partially with respect to x, yielding where Evaluation of Equation 8 at x = x 0 , with C given by Equation 7, produces the result (by choice of C), the boundary condition required to integrate Equation 4 is then clearly given by The main advantage of working with Y, as opposed to F, is that its entries never grow excessively large or small.Equation 4 can therefore be integrated without difficulty over several hundreds of nanometres, even in the presence of strong k•p coupling.
Once Y has been obtained by numerical integration of Equa- tion 4 over the entire heterostructure, the tunnelling probabilities for each channel can be obtained from the reflection matrix R, defined next.

Reflection matrix
In the left contact region (see Fig. 1), for x < x L , the solution to Equation 2 can also be expressed in terms of plane waves.In this region there are both incident and reflected waves and the solution must be written as where D± are defined as before and R is the desired reflection matrix.By evaluating Equation 12and its first derivative at x = 0, and multiplying both the resulting equations from the right by C -1 , it is found that and Here K is defined as in Equation 9. Multiplication of Equation 14 from the right by the inverse of Equation 13, produces By solving for R in Equation 15, the reflection matrix is obtained as All quantities on the right of Equation 16 are to be evaluated at x = 0.The matrix Y(0) is obtained by integrating Equation 4 from some point x = x 0 > x R back to x = 0, starting with the boundary condition given by Equation 11.Note that the diagonal entries in R contain the reflection amplitudes for each channel.It can be shown that transmission coefficient (no sum over i), where R ii is the ith diagonal entry of R and R ii * denotes complex conjugation.

Self-consistent Schrödinger-Poisson scheme
Because the above method produces the log-derivative Y(x) it can be used to perform a fully self-consistent calculation of the heterostructure energy-band profile.This is an important advantage over other methods, such as the scattering matrix method which only produces the amplitudes of the envelope functions, not the functions themselves.Once Equation 4 has been integrated and Y(x) is known throughout the entire heterostructure, a second integration can be performed to obtain the envelope function matrix F(x) from Equation 5and the boundary condition F(x 0 ) = 1.By using the envelope functions, the energy-band profile and the corresponding reflection matrix can then be calculated self-consistently in conjunction with Poisson's equation. 16

Current-voltage characteristics
Having obtained the transmission coefficients T i as functions of the electron energy E, transverse component of the wave vector k || and applied voltage V a ; the current I flowing through the device can be calculated as a function of V a .For the InAs/GaSb/InAs quantum well device studied in Section 3, there is only one open channel and the current I flowing through the structure is defined as 17,18 where FA and FC are the respective Fermi distribution functions for the left (anode) and right (cathode) contact regions.In Equation 17, A is the device area, T is the transmission coefficient, ρ is the density of states and νx is the incident carrier velocity.The latter two quantities are computed numerically by using the appropriate bulk carrier energy-dispersion relation.

Results
0][21] In this model the basis has been chosen in order to decouple the spin components of the electron for k z = 0.Each of the two possible spin states of the electron can then be treated separately in terms of two identical 3×3 matrix Hamiltonians.For details concerning this model, including the k•p parameter values for InAs and GaSb, see Botha. 4 There are a variety of numerical methods available for integrating coupled systems of first order differential equations such as Equation 4.Although the log-derivative form of Equation 4 is much simpler to integrate than the envelope function appearing in Equation 1, the integration subroutine nevertheless has to be quite sophisticated to take into account the discontinuities in the material parameters and potential as the integration proceeds through the heterostructure. 22The computations in this Section make use of an excellent integration subroutine called ODE.This subroutine is freely available and is described fully by Shampine and Gordan. 23he current I as a function of applied voltage V a has been computed for an InAs/GaSb/InAs quantum well device in which the GaSb layer is 10 nm wide.This device has been manufactured and studied experimentally. 3Based on the description of the device, the potential profile was calculated according to the self-consistent method outlined in Section 2.4. Figure 2 shows the calculated potential profile under zero bias and for V a = 0.6 V, respectively.In order to compute I, as defined in Equation 17, the calculated potential profile was used to obtain the corresponding transmission coefficient T according to the method described at the end of Section 2.3.Note that T is in general a function of the electron energy E, the perpendicular component of the wave vector k y (a good quantum number) and the applied voltage V a .As an example, the transmission spectrum is shown in Fig. 3 for k y = 0 and various applied voltage.The transmission spectrum differs considerably from that predicted by a simple two-band model in which the potential varies linearly. 24Instead of one broad resonance peak, corresponding to a confined light-hole resonant state in the GaSb layer, two narrower peaks are observed in Fig. 3. Physically, this more complicated transmission spectrum is the result of the coupling between a confined light-hole state in GaSb and quasi-bound electron states in the InAs anode.For k y ≠ 0 the transmission spectrum is further complicated by coupling of heavy-hole confined states in GaSb with quasibound electron states in InAs.
It is important to note that, in Fig. 2, for zero applied voltage there are two triangular quantum wells formed in the InAs layers which surround the GsSb barrier.As the applied voltage increases the potential profile gradually changes so that the triangular quantum well to the right of the barrier eventually disappears completely, as shown in Fig. 2 for V a = 0.6 V. On the other hand, an increase in applied voltage tends to increase the depth of the left triangular quantum well (in the anode).Coupling between the quasi-bound electron states in the InAs anode and the confined hole states in GaSb is responsible for the observed negative differential resistance.This coupling also accounts for the reduction in transmission peak height with increasing applied voltage, as seen in Fig. 3.
The transmission spectrum, calculated as in Fig. 3, produces the correct current-voltage behaviour for the device.Figure 4 shows the calculated room-temperature current through the device as a function of applied voltage.In this figure the calculated peak height has been normalised (using a multiplicative factor of 1.23) to the experimental peak height in the forward biased direction.This normalisation factor is required mainly because of the uncertainty in the device area (which is in the range of 10 -5 -10 -4 cm 2 ).In addition, there may be inevitable impurities in the sample and the value of the valence band offset between the two materials, and the position of the Fermi levels may also be uncertain.
The discrepancies between the current predicted by the present model and the much simpler two-band model, which was originally used to interpret the experimental data, do not alter the main conclusion which may be drawn from this experiment, namely that a resonant interband coupling mechanism is responsible for the observed negative differential resistance. 3owever, unlike the two-band current-voltage model, the present model provides semi-quantitative predictions of the peak-tovalley ratios.Since the present model does not take into account inelastic transport mechanisms, it may be concluded from the agreement reached in Fig. 3, that contrary to what was previously thought, inelastic transport mechanisms do not contribute significantly to the observed valley currents in this device, even at room temperature.Whether this conclusion generalises to other similar structures depends on how well the model compares to the other available experimental data (for example see Collins et al. 25 ).A more detailed investigation of the negative differential transport mechanism in GaSb/InAs heterostructures is currently in progress.
Finally, it should be emphasised that the preliminary results presented in this section are merely an example of how the general theory may be implemented and compared directly to an experiment.In this article, the main aim was to communicate the very general method that was developed in Sections 2.1-2.4.Since this method was based on the well-established multiband k•p model, it will be useful for theoretical studies of a wide variety of other interesting phenomena which occur in semiconductor heterostructures, for example, anisotropy (band warping) 26 and spin-splitting (bulk, structural and due to an applied magnetic field). 27

Conclusion
An alternative method was developed and implemented numerically to calculate electronic transport properties of type-II heterostructures.As an example, a realistic numerical implementation of the method was made by using a six-band matrix Hamiltonian in which the electron spin components decoupled.The preliminary results from this model have been compared to measured current-voltage characteristics of an InAs/GaSb/InAs quantum well device in which the GaSb layer width is 10 nm.The calculated current as a function of applied voltage was found to be in semi-quantitative agreement with the   Research Letters experiment 3 and it was therefore concluded that inelastic transport mechanisms probably did not play a significant role in the device.
Given that only the bulk k•p parameters for the relatively simple six-band model were used in these calculations, the semi-quantitative agreement obtained from this model is encouraging and warrants further investigation.Such an investigation is in progress.The minor discrepancies between the model and experiment considered, in particular the necessity of using a current normalisation factor, can be attributed mainly to the uncertainty in the device area and the value of the valence band offset parameter.
general relation (which holds for any non-singular matrix F) into the left-hand side of Equaion 2, the multiband k•p Riccati equation is obtained.Here Y has been defined as In recognition of the logarithmic derivative form of Equation

Fig. 1 .
Fig. 1.Schematic multiband potential profile (showing only two bands) along the growth direction of a semiconductor heterostructure.Within the active region, i.e. between X L and X R , the potential may vary arbitrarily.To the left and right of the active region there are two flat-band contact regions.V a is the applied (bias) voltage.

Fig. 2 .
Fig. 2. Calculated confining potential profile (active region only) of an InAs/GaSb/ InAs quantum well device in which the GaSb layer is 10 nm wide.The solid and dashed lines correspond to applied voltages of 0.0 and 0.6 V, respectively.

Fig. 3 .
Fig. 3. Calculated transmission coefficient T as a function of electron energy E for the InAs/GaSb/InAs quantum well device.Each of the four curves represents a particular value of applied voltage, as indicated by the legend.

Fig. 4 .
Fig. 4. Comparison of the calculated and experimental current-voltage characteristics for the InAs/GaSb/InAs device.Because the device area is uncertain, the current shown here has been normalised to the peak current in the forward bias direction.A multiplicative normalisation factor of 1.23 was used for a nominal device area of 2 × 10 -5 cm 2 .