Efficient method for transport simulations in quantum cascade lasers

An efficient method for simulating quantum transport in quantum cascade lasers is presented. The calculations are performed within a simple approximation inspired by Büttiker probes and based on a finite model for semiconductor superlattices. The formalism of non-equilibrium Green's functions is applied to determine the selected transport parameters in a typical structure of a terahertz laser. Results were compared with those obtained for a infinite model as well as other methods described in literature.


Introduction
Quantum cascade lasers (QCLs) are relatively new semiconductor devices [1]- [3].Nevertheless, they have been successfully used both in scientific research [4][5] and industry applications [6] due to their properties such as high optical power output, wide tuning range of the emitted radiation and room temperature operation.QCLs have proven to be very useful for spectroscopic applications such as remote detection of dangerous gases [7], or pollutants in the atmosphere [8].QCLs have been also proven effective in medical diagnostics where they are used as breath analyzers [9] or to study plasma chemistry [10].Scientific studies, both theoretical and experimental, have contributed much to extending the range of their use.Many different approaches can be taken to investigate QCL theoretically.For example, the Monte Carlo method was applied in [11] to predict the current density, population inversion, and optical gain.The method used by other authors [12,15] to calculate quantitatively QCLs parameters was a non-equilibrium Green's function method (NEGF) operating in a real space.Another known method based on Wannier functions is applied to QCL simulations [16,17].Infinite geometrical dimensions of the superlattice structure are assumed and periodic delocalised Bloch functions over the whole structure are considered there.It should be noted, however, that the Bloch theorem was created to describe phenomena occurring in crystal lattices, where the number of cells is incomparable with the number of semiconductor layers in the superlattice heterostructurem which are usually of several dozen, so it is difficult to properly assess the extent to which the assumed infinite dimensions of the model may affect the accuracy of QCL simulations.In [18] the authors compared the results calculated under their model to the results obtained with the finite geometrical dimensions model, and proposed Wannier-function approximation to be used for energies close to the lowest miniband, to demonstrate accuracy of their approach.We showed [19] differences in simulation results obtained with both superlattice models, i.e. with finite and infinite geometrical dimensions, which we called finite (FM) and infinite (IM) model respectively, whereas in [20,21] we presented a way for constructing the base of quantum states inspired by the Wannier functions base.The base we proposed used the wave functions that were the solution of the Shrödinger equation for the superlattice model of finite geometrical dimensions.
In the present study an efficient method for simulating quantum transport in QCL, based on the finite model is demonstrated.The performance of this method consists of several elements.(i) A quantum states base, inspired by Wannier functions for superlattice properties, is used which allows to limit the Hamiltonian representation solely to the space of energy [22].For such cases nano-device Hamiltonian matrices are small, so computations run fast.(ii) The scattering process is treated as just a contact described by diagonal matrix with phenomenological parameters η related to the scattering time [23].(iii) In the calculation scheme where nonequilibrium Green's function equations are solved, a simple approximation inspired by the Büttiker probes common in mesoscopic physics is applied [23][24][25][26][27][28][29].

Superlattice Model
As mentioned in the previous chapter the finite model for simulating quantum transport in QCL is used.The main elements of this model are shown in Fig. 1.It is dedicated to a quantum cascade laser, manufactured on GaAs/Al 0.15 Ga 0.85 As superlattice base [30].The illustrated case represents thermodynamic equilibrium.
The model comprises the finite number n of d length modules.Each module is composed of the layers in the arrangement 7.8 / 2.4 / 6.4 / 3.8 / 14.8 / 2.4 / 9.4 / 5.Such conducting layers produce a periodic potential, which models the energy scheme of the conduction band -E C with the difference in values between the quantum well bottom and the energetic barrier height ΔE C =150 meV Fig. 1.Schematic model for a cascade laser manufactured on GaAs/Al 0.15 Ga 0.85 As superlattice base [30].The illustrated case represents thermodynamic equilibrium.
Modelling technique for QCL is based on solving the Schrödinger equation for a variable effective mass, namely: where V(z) is the potential energy of the electron in the periodic potential resulting from the superlattice semiconductor layers.The solution to equation (1) for the z j <z<z j+1 region of constant potential and constant material composition can be written as wave functions in the form: The functions in the j th layer of the superlattice take on the form of plane counter travelling waves with , m j -electron effective mass, Etotal energy of an electron, V j -the potential energy.
The wave functions for the finite superlattice model defined by expression (2) are neither localized, nor periodic [19,20].They are grouped within the allowed energy minibands, which are determined by non-zero transmission coefficient value: where L stands for the total number of layers in the structure.For expression (3) the assumption that the electron beam concerns only the source (B L+1 = 0) holds and it determines the probability for an electron transition throughout the analyzed structure.The transmission coefficient Tr is equal to 1 for a certain discrete energy value of L-1 number.For the remaining values of energy within the minibands T ranges within 0 < Tr < 1. Physically it means that the electrons of the same miniband may have different transition probabilities for the structure, and it is believed to have a significant impact on the assessment of its transport properties.Figure 2.A shows the calculated Tr coefficient values for the wide energy range 0-0.08 eV.Here we see a very narrow energy ranges (minibands) for which this factor has a value close to 1.A more accurate picture of the results we can be seen in Figure 2.B where narrow energy range (0.0620-0.0625 eV) is showed.The computation method proposed in this paper uses the base of quantum states inspired by the Wannier functions.In our case the states are maximally localised and defined for a finite model of a nonsymmetrical semiconductor superlattice.These functions were sought according to the formula: where C is a normalization factor, which shapes the square modulus of the function V ν to the unity of the whole length of the structure, whereas ν indexes the minibands.

Fig. 3. Maximally localised quantum states representing 5 lowest minibands
It should be noted that Ψ ν (E,z) functions are not uniquely determined, as their particular embodiment is generally dependent on two complex constants, namely A 1 -the amplitude of the wave travelling from the source and B L+1 -amplitude of the wave travelling from the drain.Therefore the search for the maximally localized quantum states in the form of ( 4) is more difficult than the one where Wannier functions are concerned.

The Hamiltonian of the device
Hamiltonian matrix of the device can be written as ) for the cases of unpolarized and polarized superlattices respectively, while the matrix R H contains all effects associated with the scattering of electrons.The matrix of SL H can be written as [22]: where n is the number of the relevant period superlattice of length d, k is the wave vector, The Hamiltonian matrix of the polarized structure [ H takes the form [22]: where e is the electron charge, ξ is the value of the electric field, while variable P Q l R is calculated according to the formula: This relationship takes into account the interaction between both allowed minibands (indices μ, ν) and periods (index l) of the superlattice.The Hamiltonian matrix in our calculations included three periods of the superlattice and five maximally localized quantum states in each of them.
The transformation, described in the previous chapter, from the base of not localized wave functions to the base of maximally localized quantum states allows to reduce the Hamiltonian matrix representation from the energy-position representation to energy representation.This process is also advantageous in reducing the Hamiltonian matrix size (here 15x15), which improves the efficiency of the presented method.
In NEGF theory scattering processes expressed in R H are applied in the form of self-energies as described in the next section.

Green's functions and self energy approach
For calculating transport characteristics in QCL it is necessary to use nonequilibrium, i.e. retarded G R and correlation G < , Green's functions to define the nonequilibrium stationary state of the system.NEGF theory contains a description of their use for that purpose, and a self-consistent Born approximation [16] remains to be a quite popular calculation method.In this paper we propose an approach based on other methods to determine Green's functions.We were inspired by Büttiker probes (see Fig. 4), widely used in mesoscopic physics, where the parameter 6 6 * b b b i [27] describing the exchange of electrons between the quantum state of each module b and every reservoir of electrons in active device, taking into account the considered scattering effects, are introduced.These effects are included in the self-energy matrix 6 .For such an approach the matrix of the spectral function can be written in the form [27]: where N indexes all electron reservoirs.It should be mentioned here that diagonal elements of the matrix A give information about the local density of states defined as [27]: The transmission coefficient of electrons between the two reservoirs of m and n is determined with the following formula [27]: With the known local density of states, and scattering effects taken into account by the factor Γ, we can define the local concentration of electrons as [27]: where F is the function of distribution of electrons, determined by the Fermi levels μ N calculated for Büttiker probes and normalized by a factor T k B .Similarly, the current flowing through the probe to the quantum state b may be expressed as [27]: phenomenological parameter related to the scattering time ( W K / ! ) [22].The next step is to solve two kinetic equations.The first of them is the Dyson equation [23]: where R G is the matrix of the retarded Green's functions, which can be define as [24]: The second is the Keldysh equation [31] which allows to determine the correlation function It is worth emphasizing that in the method described here, the kinetic equations ( 17) and ( 19) need to be solved only once, unlike in the Born approximation method, where they are solved repeatedly in the selfconsistent process together with the equations for the self energies.
According to the assumptions made here, the density of current flowing through the laser structure, can be written just as it was for the description of the Hamiltonian as J = J 0 + J R [16], where J 0 is the current with the effects of electron scattering excluded, while J R is a component taking then into account.As shown by the authors in [32] for the structure under consideration the J R component disappears, and thus the total current density can be approximated by the J 0 defined as [32]: where the matrix > @ z H , ˆ0 is formed in the base of maximally localized quantum states variable ϑ is the volume of polarized structure.By this relationship the density of current were calculated and the obtained results are described in the next section.

Calculations of Green's functions
The model described in the previous chapters produced the simulations results illustrated below.In Figure 5 (left side) the imaginary diagonal elements of the matrixes G R and G < are plotted for wave vector k=0.There, the respective density of states (N 1D (E)) and the occupation of these states by electrons (n 1D (E)) are illustrated.The simulations were carried out with the eclectic field of 50 mV/period applied.The Fig. 5 (right side) contains also, the results of the Green's functions calculations, using the infinite model of superlattice [33].This juxtaposition shows differences in distribution of peak amplitudes of Green's functions.This is a consequence of different boundary conditions used in both models.By summing up one dimensional density of states functions for all states over the selected range for the wave vector k, the function of the density of states for the three geometric dimensions (3DOS) is obtained.This function can be written as . Fig. 6.Total density of states per period N 3D (E) calculated with FM, and occupation functions per period n 3D (E) calculated with Büttiker probes, and the eclectic field of 50 mV/period applied.The simulations were carried out for the finite and the infinite models of superlattice.simulationswere carried out for the finite and the infinite models of superlattice.In a similar way we obtained the occupation function (n 3DOS (E)) for the three geometric dimensions 6 shows graphs of N 3DOS (E) and n 3DOS (E) functions with the electric field of 50 mV/period applied.There are the results for the finite and the infinite models of superlattice.
The maps shown in Figure 7 are an extension of the results p resented in Figures 5 and 6.The density of states (Fig. 7A) and the occupation functions (Fig. 7B), defined as respectively, are illustrated.
All these functions are shown for Büttiker probes approach, with the bias voltage of 50 mV/period applied.
The results are plotted for the finite and the infinite models of superlattice, in a way that allows them to compare.The method presented is advantageous in terms of the effectiveness of simulating U-I characteristic for QCLs.An example illustrated in Fig. 8, where the current density Jo for FM and IM models was calculated according to (20) and versus voltage polarity per period is plotted.The calculations were conducted for a superlattice temperature of 25 K [30].The obtained characteristics show similarity to the results presented in [32], which were calculated using the Born approximation and in [30] where we can find the measured and simulated (with Monte-Carlo (MC) method) current-voltage characteristics for the structure considered here.It seems to confirm applicability of Büttiker probes for simulating terahertz lasers with the superlattice finite model.

Summary
The theory of Büttiker probes applied to efficient calculations of the transport parameters for quantum cascade lasers, based on the finite model of superlattice, was described.
The presented simulation proved the proposed approach to have given results in line with other applicable methods described in the literature.Superiority of the developed procedure for finding laser parameters as being much faster and easier to operate in comparison to the other methods reported elsewhere was demonstrated.It gives ground for further research into the proposed approach.

Fig. 2 .
Fig. 2. The calculated values of Tr coefficient for the wide of energy range (0-0.08 eV) A) and the same values for a selected narrow energy range (0.0620-0.0625 eV).

1 Ta
represent the off-diagonal couplings between the quantum levels in different periods, whereas parameters † n, are creation and annihilation states operators, respectively.

Fig. 4 .
Fig. 4. The profile of a generic mode is illustrated along with the placement of Büttiker probes.
the Fermi levels μ n and μ m of the Büttiker probes are determined from the condition of current continuity: the chemical potential at each iteration of each of the Büttiker probes in accordance with the relationship: no change of μ sB occurs.Having the chemical potential on the Büttiker probes, we can specify self energies Σ R and Σ < as diagonal matrices with elements in the form iη and iηf Bt (E), respectively, where

Fig. 5 .
Fig. 5. Examples of the density of states functions N 1D (E) and occupation functions n 1D (E), for k=0 calculated by use Büttiker probes, and the eclectic field of 50 mV/period applied.The simulations were carried out for the finite and the infinite models of superlattice.

Fig. 7 .
Fig.7.Density of states N(E,k||) (part A) and occupation functions n(E,k||) (part B), calculated for the bias voltage of 50 mV/period with Büttiker probes by use finite and infinite models of supperlattice.

Fig. 4 .
Fig. 4. Current-voltage characteristics of the QCL structure [30] obtained for different approaches.The blue solid line plotted the measured current-voltage characteristics of the studied structure..