AMMCR: Ab-initio model for mobility and conductivity calculation by using Rode Algorithm

TL;DR

This paper introduces an ab-initio module that calculates mobility and conductivity in semiconductors using Rode's algorithm, leveraging DFT inputs, and shows improved accuracy over traditional methods like RTA, with validation on CdS.

Contribution

The paper presents a new module integrating Rode's algorithm with DFT data for accurate mobility and conductivity calculations, interfaced with VASP.

Findings

Good agreement with experimental data for CdS

Outperforms Relaxation Time Approximation (RTA)

Provides a practical tool for semiconductor property prediction

Abstract

We present a module to calculate the mobility and conductivity of semi-conducting materials using Rode's algorithm. This module uses a variety of electronic structure inputs derived from the Density Functional Theory (DFT). We have demonstrated good agreement with experimental results for the case of Cadmium Sulfide (CdS). We also provide a comparison with the widely used method, the so called Relaxation Time Approximation (RTA) and demonstrated the improvisation of the results compared to RTA. The present version of the module is interfaced with Vienna ab-initio simulation package (VASP).

Equations

$$\frac{\partial\textit{f}}{\partial t}+\textbf{v}\cdot\nabla_{r}f+\frac{e\textbf{E}}{\hbar}\cdot\nabla_{k}f=\frac{\partial\textit{f}}{\partial t}\Bigr{|}_{\begin{subarray}{c}coll\end{subarray}}$$

$$\frac{e\text{E}}{\hslash }\cdot {\nabla }_{k}f=\int d{k}^{\prime }\{s(k^{\prime },k)f(k)(1-f(k))-s(k,k^{\prime })f(k)(1-f(k^{\prime }))\}$$

$$f(k)=f_0[ϵ(k)]+g(k)cos\theta ,$$

$$g_{k,i+1}=\frac{S_i(g_k,i)-\frac{eE}{\hslash }(\frac{\partial f}{\partial k})}{S_o(k)+\frac{1}{{\tau }_{el}(k)}}.$$

$$S_o(k)=\int [s(k,k^{\prime })(1-f^{\prime })+s(k^{\prime },k)f^{\prime }]d{k}^{\prime }$$

$$S_i(g_k,i)=\int X{g}_{k^{\prime },i}[s(k^{\prime },k)(1-f)+s(k,k^{\prime })f]d{k}^{\prime }$$

$$\frac{1}{{\tau }_{el}(k)}=\frac{1}{{\tau }_{ii}(k)}+\frac{1}{{\tau }_{ac}(k)}+\frac{1}{{\tau }_{pz}(k)}+\frac{1}{{\tau }_{dis}(k)}+\frac{1}{{\tau }_{alloy}(k)}+\frac{1}{{\tau }_{iv}(k)}+\frac{1}{{\tau }_{ni}(k)}$$

$$\mu =\frac{1}{3E}\frac{\int v(ϵ)D_s(ϵ)g(ϵ)dϵ}{\int D_s(ϵ)f(ϵ)dϵ},$$

$$v(k)=\frac{1}{\hslash }\frac{\partial ϵ}{\partial k}.$$

$$\sigma =ne{\mu }_e,$$

$$\frac{1}{{\tau }_{ii}(k)}=\frac{e^{4}N}{8\pi {ϵ}_0^2{\hslash }^2{k}^{2}v(k)}[D(k)ln(1+\frac{4k^2}{{\beta }^2})-B(k)],$$

$${\beta }^2=\frac{e^2}{{ϵ}_0{k}_{B}T}\int D_s(ϵ)f(1-f)dϵ,$$

$$N=N_A+N_D$$

$$D(k)=1+\frac{2{\beta }^2{c}^2}{k^2}+\frac{3{\beta }^4{c}^4}{4k^4}$$

$$B(k)=\frac{4k^2/{\beta }^2}{1+4k^2/{\beta }^2}+8\frac{{\beta }^2+2k^2}{{\beta }^2+4k^2}{c}^2+\frac{3{\beta }^4+6{\beta }^2{k}^2-8k^4}{({\beta }^2+4k^2)k^2}{c}^4$$

$$S_o=(N_{po}+1-{\text{f}}^{-}){\lambda }_0^{-}+(N_{po}+{\text{f}}^{+}){\lambda }_0^{+}$$

$${\lambda }_o^{+}={\beta }^{+}[(A^{+}{)}^{2}ln\mid \frac{k^{+}+k}{k^{+}-k}\mid -A^{+}c{c}^{+}-a{a}^{+}c{c}^{+}]$$

$${\beta }^{+}=\frac{e^2{\omega }_{po}{k}^{+}}{4\pi \hslash kv(k^{+})}(\frac{1}{{ϵ}_{\infty }}-\frac{1}{{ϵ}_0})$$

$$A^{+}=a{a}^{+}+\frac{(k^{+}{)}^2+k^2}{2k^{+}k}c{c}^{+},$$

$$N_{po}=\frac{1}{exp(\hslash {\omega }_{po}/k_{B}T)-1}.$$

$$S_i=(N_{po}+\text{f}){\lambda }_i^{-}{g}^{-}+(N_{po}+1-\text{f}){\lambda }_i^{+}{g}^{+}$$

$${\lambda }_i^{+}(k)={\beta }^{+}[\frac{(k^{+}{)}^2+k^2}{2k^{+}k}(A^{+}{)}^{2}ln\mid \frac{k^{+}+k}{k^{+}-k}\mid -(A^{+}{)}^2-\frac{c^2(c^{+}{)}^2}{3}]$$

$$\frac{1}{{\tau }_{ac}(k)}=\frac{e^2{k}_{B}T{E}_D^2{k}^2}{3\pi {\hslash }^2{c}_{el}v(k)}[3-8c^2(k)+6c^4(k)],$$

$$\frac{1}{{\tau }_{pz}(k)}=\frac{e^2{k}_{B}T{P}^2}{6\pi {ϵ}_0{\hslash }^{2}v(k)}[3-6c^2(k)+4c^4(k)]$$

$$\frac{1}{{\tau }_{dis}(k)}=\frac{N_{dis}{e}^{4}k}{{\hslash }^2{ϵ}_0^2{c}_l^{2}v(k)}\frac{1}{(1+\frac{4k^2}{{\beta }^2}{)}^{3/2}{\beta }^4}$$

$$\frac{1}{{\tau }_{alloy}(k)}=\frac{3\pi k^2}{16{\hslash }^{2}v(k)}{V}_0{U}_{alloy}^2\chi (1-\chi )$$

$$\frac{1}{{\tau }_{iv}(k)}=(N_e+1-f^{-}){\lambda }_e^{-}+(N_e+f^{+}){\lambda }_e^{+}$$

$$N_e=\frac{1}{exp(\frac{\hslash {\omega }_e}{k_{B}T})-1}$$

$${\lambda }_e^{+}=\frac{e^2{D}_e^2(Z-1)k{k}^{+}}{2\pi \rho \hslash v(k){\omega }_e}$$

$${\lambda }_e^{-}=\frac{e^2{D}_e^{2}Zk{k}^{-}}{2\pi \rho \hslash v(k){\omega }_e}$$

Full text

††journal: Computer Physics Communications

\cprotect

PROGRAM SUMMARY

Manuscript Title: AMMCR: Ab-initio model for mobility and conductivity calculation by using Rode Algorithm

Authors: Anup Kumar Mandia, Bhaskaran Muralidharan, Jung-Hae Choi, Seung-Cheol Lee and Satadeep Bhattacharjee

Program Title: AMMCR

Program obtainable from: www.ikst.res.in/ammcr

Journal Reference:

Catalogue identifier:

Licensing provisions: Open source BSD License

Programming language: C++

Computer: Any computer with GNU GCC compiler installed. Program is tested with gcc version 5.4.0.

Operating system: Unix/Linux

Keywords: Ab-initio, Mobility, Conductivity, DFT.

Nature of problem: Long-range interactions between electrons and longitudinal optical (LO) phonons are a major challenge in computing transport in polar semiconductors. Due to such LO phonons or generally polar optical phonons (POP), electron phonon scattering requires special attention since they can not be studied using a simple frame work such as relaxation time approximation (RTA).

Solution method: We have developed a code that calculates mobility and conductivity by using abinitio inputs with Rode algorithm which treats the interaction between the polar optical phonons and electrons in a proper way.

1 Introduction

With the ongoing development in the field of density functional theory (DFT), tremendous progress has been made in obtaining details in the electronic structure of complex materials important for understanding the physical properties of the materials in the ground state as well as in the excited state. This resulted in many computational tools which rely on first-principles approach based on DFT. To name a few: Vienna Ab initio Simulation Package (VASP) [1]-[5], Quantum Espresso [6], Spanish Initiative for Electronic Simulations with Thousands of Atoms (SIESTA) [7, 8] etc. The use of DFT obtained results as inputs for the study of the non-equilibrium properties such as transport is another area which is developing rapidly for a decade or so. Methods are being developed to study the transport properties at various levels: From semi-classical Boltzmann level to Green-Kubo formalism, Landauer formalism, non-equilibrium Green’s function formalism etc [9, 10, 11, 12]. There are few semi classical method based on Boltzmann Transport equation (BTE) [13]-[24] for the calculation of transport coefficients which use inputs obtained from either first-principles based calculations or experiments. For example, the code BoltzTrap [25] was developed in a way that it can accept the abinitio inputs from the DFT-based calculations and is widely used in scientific community. BoltzTrap is successful in studying the transport behaviour in variety of materials ranging from conductors, inter-metallics to thermoelectrics. One limitation in such formalism is the use of so called constant Relaxation time approximation (c-RTA). It is well known that Relaxation time approximation (RTA) is not appropriate for inelastic scattering mechanisms [26]. The use of c-RTA even more simplifies the problem by assuming a universal constant relaxation time for all elastic and inelastic scattering processes.

In the recent years, there has been an enhanced focus on the next-generation semiconductors made from two or more elements, particularly one that is strongly doped or highly degenerated. In specific, the so-called III-V materials are of present attention due to their very high carrier mobility compared to silicon. One important point about these compound semiconductors is that, unlike silicon, the carrier scattering due to the polar optical phonons is very essential, especially the temperature at which the device would be operational. These scattering processes are inelastic and one therefore requires a model for carrier transport that can take proper account of such inelastic scattering from a computational point of perspective.

In the present work, we have developed a code which is based on the Rode iterative method [13, 14, 15] and uses the ab-initio inputs and at present interfaced with the $\mathrm{VASP}$ code. Using such approach in recent studies, we have demonstrated the effect of inelastic scattering on the transport properties of n-type ZnSe [16] and AlGaAs2  [17]. The code uses electronic band structure, density of states, phonon frequencies, dielectric, piezoelectric and elastic tensors as inputs for the simulation. The band structure and density of states must be in $\mathrm{VASP}$ output format, while the other inputs such as the dielectric, elastic and piezoelectric constants could be either from abinitio calculations or from the experimental results depending on the choice of the user.

In this manuscript, we present a general implementation of AMMCR module which can be easily interfaced with VASP tool. The module AMMCR is written in C++. In this module eight types of scattering mechanisms are included, these are ionized impurity scattering, polar optical phonon scattering due to longitudinal phonon, acoustic deformation scattering, piezoelectric scattering, dislocation scattering, alloy scattering, intra-valley scattering and neutral impurity scattering. Out of these eight scattering mechanism any scattering can be included or excluded from simulation. We have tested AMMCR for CdS and found reasonable agreement with the experiment. We have also demonstrated that results with the Rode iterative method is much better than relaxation time approximation (RTA).

This paper is organized as follows: section 2 describes the the Boltzmann transport formalism, Rode algorithm and different scattering mechanisms; section 3 describes the structure of the code; section 4 describes the input files required for simulation; test system for Cadmium Sulfide (CdS) is discussed in section 5 and we conclude in section 6.

2 Theoretical framework

2.1 Solution of Boltzmann Transport Equation via iterative method

The motion of the carriers are described via a probability distribution function f(r,k,t) in real and momentum space and as a function of time. f(r,k,t) is the solution of the Boltzmann transport equation given by [26, 27, 28]

[TABLE]

Where v is the group velocity of the carriers, E is the applied electric field, \frac{\partial\textit{f}}{\partial t}\Bigr{|}_{\begin{subarray}{c}coll\end{subarray}} represents the change in the distribution function with time due to collisions.

For steady state, $\frac{\partial\textit{f}}{\partial t}=0$, and assuming spatial homogeneity ($\nabla_{r}f=0)$, the Eq.1 can be written as

[TABLE]

$S(k^{\prime},k)$ represents scattering rate for an electron making a transition from a state $k$ to another state $k^{\prime}$. In the presence of low electric field the distribution function is written as

[TABLE]

Where $f_{0}[\epsilon(k)]$ is the equilibrium distribution function, and $\cos\theta$ is the angle between applied electric field and $k$. $g(k)$ is the perturbation to the equilibrium distribution function. Higher order terms are neglected here, since transport coefficients are calculated under low electric field conditions. Now, we have to calculate the perturbation $g(k)$ to calculate transport coefficients. The perturbation in distribution function can be obtained by inserting the $f(k)$ and $f(k^{\prime})$ from the Eq.3 to the Eq.2. After some manipulation, $g(k)$ can found as [13, 14, 15, 16, 17]

[TABLE]

As $g_{k}$ appears on both the left and right hand sides of the Eq.4, the above equation should be solved in iterative manner. Here i is the iteration index. $\frac{1}{\tau_{el}(k)}$ is the sum of momentum relaxation rates of all elastic scattering rates, $S_{i}$ and $S_{o}$ are in scattering and out scattering operator for inelastic scattering and are given by [15].

[TABLE]

[TABLE]

where X is the cosine of the angle between the initial and the final wave vectors and $f^{\prime}=f(k^{\prime})$. Since convergence is exponential [15], it requires only a few iterations for convergence. $\frac{1}{\tau_{el}(k)}$ is given by

[TABLE]

where the subscripts el,ii, ac, pz, dis, alloy, iv and ni are stands for elastic, ionized impurity, acoustic deformation potential, piezoelectric, dislocation, alloy, intra-valley and neutral impurity scattering processes respectively. For calculating mobility, thermal driving force, $v(\frac{\partial f}{\partial z})$ in Eq. 4 is set to zero, only electric driving force $\frac{eE}{\hbar}(\frac{\partial f}{\partial k})$ is considered. The carrier mobility $\mu$ is given by [13, 14, 15, 29]

[TABLE]

where $D_{S}(\epsilon)$ represents density of states. The carrier velocity $v(k)$ is calculated from abinito band structure by using relation

[TABLE]

After calculating mobility electrical conductivity can be calculated by

[TABLE]

where $n$ is the electron carrier concentration, and $\mu_{e}$ is the electron mobility.

2.2 Scattering Mechanisms

Eight different types of scattering mechanisms are included in this code. Out of these eight scattering mechanisms any scattering mechanism can be included in the simulation. These eight scattering mechanisms are ionized impurity scattering, Polar Optical phonon scattering due to longitudinal phonon, acoustic deformation scattering, piezoelectric scattering, dislocation scattering, alloy scattering, intra-valley scattering and neutral impurity scattering.

Interaction of electrons with ionized impurity potential results in ionized impurity scattering. Ionized impurity scattering is an important scattering mechanism at higher doping concentration and at lower temperatures. For ionized impurity scattering Brooks-Herring approach [30] is used. The momentum relaxation rate for ionized impurity scattering is given by [15, 29]

[TABLE]

where $\epsilon_{0}$ is the dielectric constant and $\beta$ is the inverse screening length given by

[TABLE]

where N is the concentration of ionized impurity and it is given by

[TABLE]

where $N_{A}$ and $N_{D}$ are the acceptor and donor concentrations respectively. The expressions for D(k) and B(k) are [15].

[TABLE]

[TABLE]

The wave function admixture $c(k)$ is the contribution of p-orbital to the wave function of the band. It is calculated by projecting the Kohn-Sham wavefunctions onto the spherical harmonics which are non-zero only within the spheres centering the ions and this is already implemented in VASP package.

Polar optical phonon (POP) scattering is inelastic and an-isotropic scattering mechanism. In most of the polar semiconductors POP scattering is dominant scattering mechanism near room temperature or in the higher temperature region. The out scattering operator is given by [15]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, a(k) gives the contribution of s-electrons to the band, and in a way similar to c(k) it is calculated from the first-principles. $k^{\pm}$ is the solution of the equation, $\varepsilon(k)\pm\hbar\omega_{po}$. Any quantity with subscript plus/minus has to be evaluated at energy corresponding to k+/k-. $\lambda_{o}$ are rate of out scattering. The subscript plus and minus denotes the scattering by the absorption or emission respectively, so it is to be evaluated at an energy $\epsilon+\hbar\omega_{po}$ for absorption and at energy $\varepsilon(k)-\hbar\omega_{po}$ for emission. If the energy of phonons is less than $\hbar\omega_{po}$, the emission of phonons is not possible and hence $\lambda^{-}_{o}$ is to be considered to be zero.

$N_{po}$ is the number of phonons and is given by [15]

[TABLE]

The scattering in operator $S_{i}$ is given by [15]

[TABLE]

Here $\lambda_{i}$ is the rate of in scattering. Again, the plus and minus subscriptions indicate absorption and emission, respectively.

[TABLE]

Coupling of electrons with non-polar acoustic phonons results in acoustic deformation potential scattering. The momentum relaxation rate for acoustic deformation potential scattering is given by [15, 29]

[TABLE]

Where cel is the spherically averaged elastic constant. $E_{D}$ is acoustic deformation potential and is given by conduction band shift (in eV) per unit strain due to acoustic waves.

Polar scattering due to acoustic phonons is called piezoelectric scattering. Piezoelectric scattering is important at low doping concentration and at low temperature in polar materials. The momentum relaxation rate for piezoelectric scattering with ab-initio parameters is given by [15, 29]

[TABLE]

where $P$ is a dimensionless piezoelectric coefficient, it is isotropic for zinc blende structure and anisotropic for wurtzite structure (see the appendix).  

The momentum relaxation rate for dislocation scattering is given by [31]

[TABLE]

In case of alloy there is one more scattering mechanism due to atomic disorder, known as alloy scattering. The momentum relaxation rate for alloy scattering is given by [32]

[TABLE]

The momentum relaxation rate for intra-valley scattering is given by [15]

[TABLE]

where $N_{e}$ is the phonon occupation number and given by

[TABLE]

$\lambda_{e}^{+}$ is given by

[TABLE]

Similarly $\lambda_{e}^{-}$ is given by

[TABLE]

where $Z$ is the number of equivalent valleys, $D_{e}$ is the intervalley deformation potential (units of electron volts per meter), $\hbar\omega_{e}$ is the phonon energy, $\rho$ is the density of the material. If $E<\hbar\omega$, $\lambda_{e}^{-}$ is considered to zero.

There is one impurity scattering due to non ionized donors called neutral impurity scattering. Neutral impurity scattering is important at higher doping concentration and low temperature. For neutral impurity scattering we have used Erginsoy model[33]. The momentum relaxation rate of neutral impurity scattering is given by

[TABLE]

where $N_{n}$ is concentration of neutral impurities in semiconductor.

2.3 Ab-initio Inputs

For calculating transport properties, band structure and density of states of semiconductor material is required. The required band structure and density of states are calculated by using the density functional theory (DFT) using a three dimensional k mesh around the conduction band minimum (CBM). Near CBM much finer meshing is required to obtain good results, since this near region to CBM will play a major role at low electric field. To obtain a k-point file in reciprocal space for DFT calculation a program k_point_generator.cpp is given with the code. The program is saved in the ‘utility’ folder. The conduction band is represented by a function of distance from the CBM, by taking the average of the energy values of the k-points that are at the same distance from CBM. We have done an analytical fitting of the band with a six degree polynomial, to get smooth curves for group velocity. The group velocity is calculated by using equation

[TABLE]

The calculated group velocity by the above equation is used to calculate required different scattering rates. All other inputs required for calculations, low and high frequency dielectric constant [34, 35], elastic constant, piezoelectric constant and polar optical phonon frequency $\omega_{po}$ [36] are calculated by using Density functional theory (DFT) by using VASP. The acoustic deformation potential is given by equation

[TABLE]

where V is the volume, $V_{0}$ is the volume under zero pressure and CBM is conduction band minima [37, 38]. So, all required inputs are calculated by using DFT. Only crystal structure of semiconductor material is given as input. So, there is no need to heavily depend on experimental data for mobility and conductivity calculation.

2.4 Simulation Flowchart

Fig. 1 shows the step to calculate mobility and conductivity using AMMCR code. First, we have to calculate band structure and all other input parameters as explained in the previous section by using first principles methods. Then analytical fitting of band structure is done, to obtain a smooth curve for group velocity. The Fermi level is calculated with smooth band structure obtained after analytical fitting by using equation,

[TABLE]

where $D_{S}(\epsilon)$ represents the density of states at energy $\epsilon$, where $\epsilon_{c}$ represents the bottom of the conduction band and $V_{0}$ represents the volume of the cell. Next, scattering rates are being calculated for different selected scattering mechanisms by using Eq.s 11, 16, 23, 24, 25, 26, 27, 32. The perturbation in distribution function $g(k)$ is calculated by using Eq.4, keeping $S_{i}(k)=0$, this particular step gives the results within RTA. To obtain the results beyond RTA, $g(k)$ in is calculated by iteration till $g(k)$ (or $g(\epsilon)$) converges. The $g(\epsilon)$ obtained in the above manner is used to calculate the transport coefficients.

3 Code Layout

Upon unzipping the tar file, the main folder ”AMMCR” contains four sub-folders: ‘src’, ‘utility’,‘example ’ and ‘manual’. ‘src’ folder contains different functions to calculate mobility and conductivity. It contains a file ‘main.cpp’, it is the main file, that calls all other functions for calculation. ‘utility’ folder contains a file ‘k_point_generator.cpp’, this file contains a program that will generate a ‘k_points_file ’ file, containing the k-points required for the DFT simulation of the band structure using VASP. ‘manual’ contains the user-guide while the folder ‘example ’ contains the example files.

4 Input Files

The execution of AMMCR module requires five files as input. Four input files EIGENVAL, OUTCAR, PROCAR and DOSCAR are to be obtained by using VASP package and one input file input.dat that contains the value of different material constants calculated by abinito method as explained in section 2.3. Sample input files are given in the ‘example’ folder. Typical inputs for a given donor concentration is given below (for the case of CdS),

32.75 37.75 48.43 61.05 75.55 99.12 123.99 154.08 154.13 171.998 194.67 264.20 299.94 346.74 411.18 # Temperature loop 6.9e15ΨΨΨΨΨ # Donner conc in cm^-3 Ψ 0ΨΨΨΨΨΨ # Acceptor conc in cm^-3 Ψ 0ΨΨΨΨΨΨ # Neutral impurity conc cm^-3 Ψ 9.97ΨΨΨ # Static dielectric const 6.1ΨΨΨΨΨΨ # High frequency dielectric constΨ 0ΨΨΨΨΨΨ # Band gap (eV) 1ΨΨΨΨΨΨ # Equivalent number of valence band valleysΨ 1ΨΨΨΨΨΨ # Equivalent number of conduction band valleys 0ΨΨΨΨΨΨ # Density in gm/cm^3 1 1 1 1 0 0 0 0ΨΨ # Different scatt mechanism included 0ΨΨΨΨΨΨ # Dislocation density cm^-2 6.44ΨΨΨΨΨ # Longitudinal Optical freq in THZΨ 12ΨΨΨΨΨΨ # Acoustic deformation potential in eV 0.141104752018947 ΨΨ # Piezo Electric coefficient (dimensional less) 7.6924e+11 ΨΨ # Longitudianl elastic constant in dyne/cm^2 1.9146e+11ΨΨΨΨ # Transverse elastic constant in dyne/cm^2 0ΨΨΨΨΨΨ # Alloy potential in eV 0ΨΨΨΨΨΨ # Volume of unit cell in nm^-3 0ΨΨΨΨΨΨ # Fraction of atom for alloy 0ΨΨΨΨΨΨ # Phonon freq for intravalley scatt in THz 0ΨΨΨΨΨΨ # Coupling constant for intravalley scattering in 10^8 eV/cm 0ΨΨΨΨΨΨ # Number of final valley for intravalley scattering 0ΨΨΨΨΨΨ # Density of state 0 read from DOSCAR; 1 for free electron density 10ΨΨΨΨΨΨ # Maximum number of iteration

5 Test System

We have calculated the transport properties of bulk CdS using the AMMCR code. We have considered the wurtzite structure of CdS (space group No.:186). The electronic structure calculations were performed within the frame-work of density functional theory (DFT) with Perdew-Burke Ernzerhof exchange correlation energy functional which is based on a generalized gradient approximation. Fig. 2 shows electronic band structure and density of states (DOS) obtained for CdS. We have calculated all required input parameters by using abinito principle. As an input for the transport calculation within Rode’s method, only band structure for one valley is needed, we have therefore performed non-self consistent calculations of the band energies in a special k-point mesh around the $\mathrm{\Gamma}$ point with 8531 k-points. Such an approach enables us to efficiently account for the group velocity. As CdS has wurtzite structure, it has different dielectric permittivity in parallel and perpendicular to the c axis. For calculation isotropic value of permittivity given by equation

[TABLE]

is used. Calculated value of low and high frequency dielectric constant are 9.97 and 6.1 respectively. We have obtained a abinitio value of polar optical phonon frequency of 6.44 THz. We have obtained abinito value of lattice constant $a_{0}=0.420nm$ and $c_{0}=0.687nm$, 12 eV acoustic deformation potential, 0.10059 and 0.141104 of perpendicular and parallel piezoelectric coefficient. Calculated value of elastic constant are $c_{11}=7.77\times 10^{10}N/m^{2}$ , $c_{13}=3.72\times 10^{10}N/m^{2}$, $c_{33}=8.79\times 10^{10}N/m^{2}$ and $c_{44}=1.49\times 10^{10}N/m^{2}$ respectively.

Fig 3 and 4 shows scattering rate as a function energy for different doping concentration $N_{D}=1\times 10^{13}cm^{-3}$ and $N_{D}=1\times 10^{15}cm^{-3}$ at temperatures 30 K and 300 K. Piezoelectric scattering is considered to be the significant scattering mechanism at lower temperature and lower doping concentration. At 30 K average energies of carriers is $\frac{3}{2}k_{B}T=0.0038eV$, so most of carriers lies in the low energy region. For CdS at lower temperature the piezoelectric scattering is most dominant scattering mechanism as in fig 3(a) for doping $N_{D}=1\times 10^{13}cm^{-3}$ while for doping $N_{D}=1\times 10^{15}cm^{-3}$ ionized impurity and piezoelectric scattering are most dominant scattering mechanisms as in fig 4(b). While at room temperature for both doping concentrations polar optical phonon scattering is the most dominant scattering mechanism for CdS as it is clear in both figures Fig. 3(b) and fig 4(b). So from room temperature to the higher temperature region POP scattering is most dominant one. In both Fig. 3 and Fig. 4 there is a sudden change in POP scattering rate after a particular energy, this is due to the fact that if an electron energy is less than POP energy $\hbar\omega_{po}=0.0266eV$, the electron can scatter only by the absorption of optical phonons and if the electron energy is higher than energy $\hbar\omega_{po}$, then it can scatter through both absorption and emission of phonons.

We have calculated mobility and conductivity of CdS using AMMCR module for different temperature and doping concentrations. Fig. 5 shows mobility and fig 6 show conductivity as a function of temperature for doping $6.9\times 10^{15}cm^{-3}$ [39]. Fig 5 shows good qualitative and quantitative agreement between experimental and calculated curves. Fig. 5 and fig 6 shows the results obtained with both Rode algorithm and RTA. Since piezoelectric constant is different in parallel and perpendicular to c-axis of CdS crystal, so there are two theoretical curves shown for both Rode and RTA. At lower temperature there is a deviation between experimental and theoretical curves, this is due to the presence of neutral impurities in the sample. Rode curve has an average relative error of 10.93 % and RTA curve has an average error of 37.06 % above 90 K with parallel piezoelectric coefficient. So, Rode results are much better than RTA results. This is due to POP scattering, since POP is inelastic as well as an-isotropic scattering mechanism, so it modulated electron energy and it is most dominant scattering mechanism around room temperature to higher temperature region, so RTA is inappropriate for it. At a lower temperature, both RTA and Rode show about the same mobility. This is due to the fact that POP scattering is weak at a lower temperature.

As already mentioned, most of the compound semiconductors are polar in nature. In a polar semiconductor, usually polar scattering is the most dominant scattering mechanism.Therefore, AMMCR will serve a better purpose in comparison to the RTA based codes such as Boltztrap. Fig. 7 shows the contribution of mobility from different scattering mechanisms. Above 70 K polar optical phonon scattering is most dominant scattering mechanism, below 55 K ionized impurity scattering is most dominant scattering mechanism and in the intermediate region acoustic deformation scattering is most dominant scattering. Fig. 8 and fig. 9 shows mobility variation with doping for temperature 77 K and 300 K respectively by assuming a compensation ratio of unity. The curves shows decrease in mobility with increasing doping concentration, this is due to increase in number of ionized centers with increasing doping. In the Fig.10, we show the convergence of the non-equilibrium part of the distribution function. It can be seen that the convergence is achieved after a few iterations.

6 Conclusion

We presented a module AMMCR to compute mobility and conductivity of semiconductor materials. The module is written in C++ and is currently interfaced with VASP. We have tested the code for CdS and obtained good agreement with the experimental results.

7 Acknowledgments

The authors acknowledge funding from Indo-Korea Science and Technology Center (IKST). We thank the Computing infrastructure from the Center of Excellence in Nanoelectronics, IIT Bombay and IKST for providing the computational facility.

8 Appendix

For zinc blende structure piezoelectric coefficient is given by [15]

[TABLE]

where $h_{14}$ is one element of piezoelectric stress tensor and $c_{l}$, $c_{t}$ are the spherically averaged elastic constant for longitudinal and transverse modes respectively and are given by [15] equations

[TABLE]

[TABLE]

where $c_{11}$, $c_{12}$ and $c_{44}$ are three independent elastic constants. For wurtzite structure we use piezoelectric coefficients $P_{\parallel}$ and $P_{\perp}$ for mobility measured with electric field parallel and perpendicular to the c axis of the crystal. For wurtzite structure piezoelectric coefficients $P_{\parallel}$ and $P_{\perp}$ are given by [15]

[TABLE]

[TABLE]

[TABLE]

where $h_{15}$, $h_{31}$ and $h_{33}$ are the three independent elements of the piezoelectric stress tensor of wurtzite structure and $c_{l}$ and $c_{t}$ are spherically averaged elastic constant, there are given by equations [15]

[TABLE]

[TABLE]

where $c_{11}$, $c_{13}$, $c_{33}$ and $c_{44}$ are elastic constants.