The numerical method for solving the transport equations in HgCdTe heterostructures using the non-equilibrium distribution functions

This paper presents a numerical method for solving the set of transport equations in semiconductor heterostructures by using the non-equilibrium distribution function for electrons and holes. This method enables the calculation of carrier concentration, carrier mobility and entropy by integrations in the space of wave vectors. In the same way the electrical current density and density of entropy currents are determined. The influence of quasi-Fermi energy gradients for electrons and holes and the gradient of temperature on the physical parameters of the heterostructure is taken into account.


Introduction
Usually semiconductor devices work in non-equilibrium conditions caused by electrical biasing and optical excitations. The flows of electric carriers and energy are driven by gradients of intense state parameters; quasi-Fermi energies for electrons and holes rU n r ð Þ and rU p r ð Þ respectively, and the gradient of temperature rT. These gradients can have a significant impact on the non-equilibrium distribution functions determining the statistical This article is part of the Topical Collection on Numerical Simulation of Optoelectronic Devices 2016.
Guest edited by Yuh-Renn Wu, Weida Hu, Slawomir Sujecki, Silvano Donati, Matthias Auf der Maur and Mohamed Swillam. probability that the electrical carrier has position r and wave vector k. However, in commercially available simulators the non-equilibrium distribution functions do not contain these gradients, while one can find other non-equilibrium approaches including standard hydrodynamic transport and hydrodynamic transport with six moments of the Boltzmann transport equation (for example Grasser et al. 2001). Typically, the nonequilibrium distribution functions for electrons and holes with these gradients are derived by solving the kinetic Boltzmann equation, as it was done for example in the classical work (Marshak and van Vliet 1984). We have shown (Jóźwikowski et al. 2015), that functions derived in such way are only an approximation of more general terms. To derive them we had postulated a modified variational principle for the semiconductor heterostructure in steady state conditions being the composition of the Massieu function with the Gyarmati's principle. This proposed principle for conduction band electrons in steady state conditions reads as: here s e is the density of electron entropy, n is electron concentration, u e is the density of electron energy, T the temperature,U n is the quasi Fermi energy for electrons, s e is the mean electron relaxation time, r s e is the efficiency of the source of electron entropy, V is the volume of the device and W e is the dissipation function in force representation (see Gyarmati 1967Gyarmati , 1969Verhas 2014;Onsager 1931a, b;Onsager and Machlup 1953 The solution of Eq. (2) is the following non-equilibrium distribution function for electrons: A similar formula holds for non-equilibrium distribution function for holes, derived by this same method as for electrons: In this paper distribution functions for ionized donors N D + and acceptors N À A are expressed in a standardized form [see for example (Blakemore 1962)]. In relations (3) v r ð Þ denotes electron affinity, m e * it's effective mass, e elementary charge, W electrical potential, T temperature, s e density of electron entropy, ℏ Planck's constant, k B Boltzmann's constant, n electron concentration, U n quasi-Fermi energy for electrons, E g is energy gap, s e denotes mean relaxation time, ris the position vector and k is the wave vector. Subscripts p and h in Eq. (4) relate to the holes. Electron kinetic energy e k e k and heavy hole kinetic energy e k h k are expressed by Kane's relations (Kane 1957). Due to the scalar products k Á rU n , k Á rU p and k Á rT the functions expressed by Eqs. (3) and (4) are not even in k. However in the analyzed structures ∇T is very small, so also expressions containing it can be omitted. Therefore integrating their product with the k in the first Brillouine's zone (BZ) allows to determine the density of carrier current and the density of entropy current. In the case of HgCdTe two different scattering mechanisms, for electrons and holes, i.e. scattering on ionized impurities as well as scattering on polarized optical phonons (Nag 1980;Lundstrom 1990) mainly limit the relaxation time. We have presented the detailed description of relaxation time for electrons and holes in HgCdTe in our previous paper (Józwikowski et al. 2004). If we consider the small values of gradients of rU n , rU p and rT and expand relations (3) and (4) into Taylor series one obtains these same expressions as those obtained by (Marshak and van Vliet 1984). The physical parameters of HgCdTe can be found at (Capper and Garland 2011).

Method of numerical analysis
Our goal is to determine the physical parameters of devices, through the appointment of non-equilibrium distribution function. To reach it one has to solve the following non-linear system of transport equations, namely Poisson equation, continuity equation for electrons and holes, and the equation of energy balance, with suitable boundary conditions on the detector's surface and its electric contacts (see for example Parrot 1996; Jóźwikowski et al. 2010). In this paper we have considered Ohmic contacts and Neumann's boundary conditions on the surface beyond contacts. We have assumed constant temperature on the detector's surface and electric contacts, the same as in the environment.
The numerical solution of the set (5, 6, 7, 8) is obtained by using the Newton iterative methods. For this purpose we use our own computer programmes where the cell-centred finite volume method (FVM) as a discretization method is applied. G opt denotes the optical generation rate,ẽ k e;opt andẽ k h;opt is the average kinetic energy generated optically of electrons and holes, respectively, c V is the specific heat,v is the thermal conductivity coefficient, G is the carrier generation rate, R is the carrier recombination rate, ɛ 0 is the permittivity of free space and ɛ is the dielectric constant. The G and R terms include interband Auger1 and Auger 7 mechanisms (Landsberg 1991). In addition, we have included the Shockley-Read-Hall processes caused by vacancies of metals and dislocations. A detailed description of the expressions for G and R being the function of quasi-Fermi energies and temperature can be found in our works (Jóźwikowski et al. 2010(Jóźwikowski et al. , 2012. G opt is omitted in this paper. Other quantities appearing in set (5,6,7,8) are calculated by numerical integration in spherical space of wave vector k (Lundstrom 1990). For electrons, these values are expressed as follows: Here j e r ð Þ denotes the density of electron current and j s e r ð Þ is the density of electron entropy current. Similar expressions for quantities related with holes take forms: We count the individual components of the current density in the standard way, for example x-th component of j e ; reads as The upper limit for the module of k vectors in our calculations for integrating in k space [in relations (9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19)] is the value of k max = 10 10 m −1 and it is due to our practical experience.
While performing numerical simulations of semiconductor devices, we assume the knowledge of the spatial distributions of molar composition, electron affinity and doping. For those distributions we use the local electrical neutrality condition, which enables us to determine the initial value of the electric potential W 0 r ð Þ. We then solve Poisson equation using iterative methods in order to determine the iterative corrections dW r ð Þ for the electric potential in thermal equilibrium (Jóźwikowska 2008). This, in turn, allows the calculation of n, p, s e and s h by using relations (9), (14), (10) and (17), respectively. In conditions of thermal equilibrium distribution functions for electrons and holes are expressed by relations (3) and (4), respectively, but in this case U n r ð Þ ¼ U p r ð Þ ¼ F, where F is the Fermi energy and there is no gradients of intense state parameters. In biased heterostructures these gradients occur. Now, to determine the change of W r ð Þ; U n r ð Þ; U p r ð Þ; T r ð Þ, with the change of the bias voltage, we have to solve the set of Eqs. (5,6,7,8). These equations are solved now by fully coupled Newton iterations. It allows us to determine the iterative corrections dW r ð Þ; dU n r ð Þ; dU p r ð Þ and dT r ð Þ: Next, the corrected parameters W r Eqs. (9,10,11,12,13,14,15,16,17,18) to calculate the corrected values of f e , n, μ e , j e , s e , j s e , f h , p, μ h , j h , s h and j s h .

Some numerical results
Below we present results of calculations of physical parameters of the cylindrical Hg 1- x Cd x Te barrier MESA detectors with a radius of 30 μm (length of A'B line on Fig. 1), working at 300 and 230 K. Spatial distributions of mole fraction and donor and acceptor concentrations are shown along the axis of symmetry of the structure (line AA' in Fig. 1) in Fig. 2. All the spatial distributions of physical parameters of the structure are applied to distributions along AA' line. In the calculations we assumed that the dislocation density in the whole structure is equal 10 6 cm −2 except for areas where there is a gradient of molar composition. In these regions the density of dislocations is increased according to the relationship proposed by (Yoshikawa 1988). The density of metal vacancies, being trap centres in Hg 1−x Cd x Te, is assumed to be constant in all the heterostructures. We have considered their concentrations equal to 10 14 cm −3 and 10 15 cm −3 . Figure 3 shows calculated spatial distribution of carrier concentration and their mobility at 300 K in structure biased with 1 V in reverse direction. Spatial distribution of calculated energy band diagram and quasi Fermi energy for electron Φ n and holes Φ p , are shown in Fig. 4. Experimental and theoretical calculations of normalized current-voltage characteristics are compared in Fig. 5a, b. The calculations were performed taking into account the existence of serial resistance. Measured, normalized value of it at 300 K is equal 0.05 Ω cm 2 . Figure 5a shows the results at 230 K. Good agreement between theoretical and experimental results is obtained for trap concentrations 10 14 cm −3 ≤ N T ≤ 10 15 cm −3 . Figure 5a shows results at 300 K. However, detailed considerations need to take into account different trap concentration in different regions of the heterostructure. This is probably the reason of discrepancies between experimental and theoretical results for bias voltage in range 100-400 mV.

Conclusion
We have applied, derived earlier by ourselves, the non-equilibrium distribution functions for electrons and holes in the numerical modelling of barrier MESA infrared detectors. These functions are more general than those commonly obtained by solving the BKE and are valid also for strong gradients of quasi Fermi energies and the gradient of temperature. The proposed numerical method enables the calculation of physical parameters by the numerical integration of distribution function in the space of wave vector k. It is especially important to obtain proper values of the current densities. The method will be developed in future to distinguish the effective temperature of electrons, holes and lattice as well as the quasi-Fermi energies for electrons occupy donor, acceptor and trap centres.