Thermoelectric transport properties of silicon: Towards an ab initio approach
We have combined the Boltzmann transport equation with an ab initio approach to compute the thermoelectric coefficients of semiconductors. Electron-phonon, ionized impurity, and electron-plasmon scattering rates have been taken into account. The electronic band structure and average intervalley deformation potentials for the electron-phonon coupling are obtained from the density functional theory. The linearized Boltzmann equation has then been solved numerically beyond the relaxation time approximation. Our approach has been applied to crystalline silicon. We present results for the mobility, Seebeck coefficient, and electronic contribution to the thermal conductivity, as a function of the carrier concentration and temperature. The calculated coefficients are in good quantitative agreement with experimental results.
Interest in thermoelectric materials has been rapidly growing in recent years. This is in part due to the new expectation of materials with a higher dimensionless figure of merit brought about by the nanotechnology revolution.Boukai:2008 (); Harman:2002 (); Hochbaum:2008 (); Venkatasubramanian:2001 () From the theoretical point of view, it is important to be able to predict various thermoelectric properties without resorting to adjustable parameters. However, despite great advances in predicting the electronic structure of materials, calculation of thermoelectric transport properties from first principles still presents a challenge, even in the case of simple bulk materials.
In some relevant works, attempts have been made to use a from first principles description of the band structure, combined with the relaxation time approximation (RTA) of the scattering mechanisms.Wang:2009 (); Oh:2008 (); Huang:2008 (); Singh:2007 (); Wilson-Short:2007 (); Hazama:2006 (); Madsen:2006 (); Madsen:2006b (); Ishii:2004 (); Thonhauser:2003 () However, these approaches can still be considered as semi-empirical, since the scattering rates depend on adjustable parameters, and moreover rely on approximations for their energy dependence. A recent calculation of the mobility of SiGe alloys has been presented in Ref. Murphy-Armando:2008, , where the scattering rates have also been computed from first principles in the RTA. Regarding silicon, its mobility has been explored employing a parameter free approach Restrepo:2009 () within the RTADziekan:2007 (); Yu:2008 (), and its lattice thermal conductivity has been computed ab initio, together with that of germanium and diamond.Broido:2007 (); Ward:2009 ()
Thus, previous ab initio works largely relied on the RTA, and in most cases, the ab initio aspect was limited to the calculation of electronic band structures. Our goal is to go beyond this stage by computing ab initio scattering rates, and solving the BTE beyond the RTA. In the present study we are focussing on the thermoelectric properties of crystalline silicon, despite the fact that it is not a good candidate for practical thermoelectric applications. However, its electronic properties have been intensively investigated in the scientific literature, with a large amount of theoretical studies devoted to the calculation of the mobility of silicon using adjustable parameters.Jacoboni:1983 (); Fischetti:1991 (); Rode1972 () Si constitutes an ideal system to test the ability of ab initio methods to predict transport coefficients.
This paper is organized as follows. The general theory for the Boltzmann equation is presented in Section II. The ab initio method used to compute the band structure and the average deformation potentials for the intervalley electron-phonon scattering is described in Section III. Results are discussed in Section IV, and conclusions are drawn in Section V.
Ii Boltzmann transport equation
The electron occupation distribution function has been obtained by solving the Boltzmann transport equation (BTE)Ashcroft:1976 (); Nag:1972 () for steady states in presence of a uniform and static external electric field and a temperature gradient ,Rode:1995 ()
where is the band energy associated with a point in the semiclassical phase space, is the chemical potential, is the electron velocity, is the unit cell volume, is the probability of a transition from to per unit of time.
The linear approximationNag:1972 () consists in assuming that can be written as a linear function of the applied external fields. At low field,
where is the equilibrium Fermi-Dirac distribution, and are the st-order corrections to due to external electric field and temperature gradient. The approximation of non-degenerate statistics is not used in this work. The combination of Eq. 1 and Eq. 2 yields two linear matrix equations,
where and , and and are -component vectors. and are matrices where is the total number of points. is diagonal with . reads , with referring to . An energy range of eV near the bottom of the conduction band has been discretized with 40 energy values. For each value of the energy , Eqs. 3 have been solved for a set of vectors on the surface of constant energy. In our calculation, discretization of the Brillouin zone is performed using the Gilat-Raubenheimer procedure (see next section).
The matrix elements of have been computed by considering different scattering mechanisms, including electronic interactions with intravalley and intervalley phonons, ionized impurities and plasmons,
where represents the electronic scattering by a phonon of vanishing wave-vector , for which, in the scattering process, the initial and final electronic states are located in the same valley. The term represents the electronic scattering by a phonon of finite wave-vector , the intervalley scattering, in which the initial and final electronic states belong to two different valleys. In silicon, this mechanism involves a scattering from the valley , located near the point, to one of the other equivalent valleys with =. The matrix represents the scattering probability of the electrons by the impurities. The latter three scattering processes , and have been calculated using models described in Jacoboni and Reggiani. Jacoboni:1983 () The deformation potentials used to compute have been calculated ab initio as described in the next section. For electronic interactions with zone-center acoustic phonons , we have used the elastic approximation, assuming the phonon energy to be negligible. The intervalley scattering is an inelastic process, with scattering probability which readsJacoboni:1983 ():
Here, is the unit cell volume, is the crystal mass density, and are the effective phonon frequency and deformation potential, characterizing a particular intervalley scattering channel, is phonon occupation number, and are the energies of the initial and final electronic states. In our calculations, the following approximation has been made:
where is the density of states. Finally, stands for the electron-plasmon interaction probability, and has been calculated using Fischetti’s approach.Fischetti:1991 ()
Once the value of the distribution function was determined, the electrical conductivity has been calculated as
The electron mobility is given by
where is the carrier concentration.
The Seebeck coefficient (thermopower) has been calculated from the solution of , using
Finally, the expression used to evaluate the electronic contribution to the thermal conductivity was
Iii Ab initio calculations
BTE is coupled with the band structure of silicon computed from first principles, and to the electron-phonon coupling constants obtained ab initio for the intervalley scattering.Sjakste:2007c (); Sjakste:2006 (); Sjakste:2007 (); Sjakste:2008 (); Tyuterev:2010 () The calculations have been performed respectively within the density functional theoryHohenberg:1964 (); Kohn:1965 () (DFT) and the density functional perturbation theoryBaroni:1987 (); Baroni:2001 () (DFPT). The local density approximation (LDA) has been used for the exchange and correlation functional. We have used the pseudopotential and plane-wave approach to solve the Kohn-Sham equations. For silicon, the pseudopotential of Refs. Giannozzi:1991, ; Tyuterev:2010, has been used, and the size of the plane-wave basis set has been limited with a cut-off energy of Ry. Finally, a Monkhorst-Pack grid has been used to sample the Brillouin zone (BZ), yielding non-equivalent points in the irreducible BZ. We have obtained an equilibrium lattice parameter of Å for silicon, and the conduction band minimum in the irreducible BZ was at , with in our calculationsTyuterev:2010 () and in the scientific literature.Madelung:1982 ()
Each isosurface has been defined by a given energy value above the conduction band edge , and has been determined from first principles. To obtain the points forming the isosurface, we have first divided the irreducible BZ into a set of small parallelepipeds, , , , where and are the unit reciprocal lattice vectors and the number of segments in the direction , respectively. In each elemental parallelepiped, a representative k-point at the given value of has been obtained with Brent’s method.Brent:1972 () The weight of the selected k-point has been computed as the area of the energy surface confined in the parallelepiped. This area has been calculated using the method proposed by Gilat and Raubenheimer.Gilat:1965 () The constant energy surface in the parallelepiped has been approximated by a plane intersecting with the parallelepiped. The plane has been chosen according to the energy gradient at the representative k-point. The details of the method can be found in Ref. Gilat:1965, . When the parallelepiped was small enough, the plane has proved to be a very good approximation to the constant energy surface.
|Jacoboni111From Table VI of Ref. Jacoboni:1983, .||This work|
|TA||0.5||140||0.61222First-order transition (see text).||138||2.88|
|LA||0.8||215||1.22222First-order transition (see text).||226||4.72|
|TO||1.26222First-order transition (see text).||698||14.56|
|LO||11.0||720||4.183330 order transition.||720||15.01|
|TA||0.3||220||0.18222First-order transition (see text).||217||4.52|
|TA||0.20222First-order transition (see text).||263||5.48|
|LA||2.0||550||1.123330 order transition.||521||10.86|
|LO||2.08222First-order transition (see text).||569||11.86|
|TO||2.0||685||2.40222First-order transition (see text).||658||13.72|
|TO||4.243330 order transition.||669||13.95|
In BTE, the implementation of the electronic scattering by the short-wavelength phonons, , has been performed with average deformation potentials computed ab initio (Table 1). Details about the -TA (transverse acoustical) and -LA (longitudinal acoustical) processes have been reported elsewhere.Tyuterev:2010 () The latter processes are forbidden by the symmetry selection rules to zero order, when the phonon wavevector is strictly equal to the vector connecting the bottoms of the two valleys or . The scattering processes are however allowed for phonon vectors that connect points in the neighbourhood of or (first-order processes). The matrix elements for these first-order scattering processes have been computed on a grid with a step of , and have been averaged by fixing the initial electronic states at the point and by summing the matrix elements over final electronic states close to the point. The above-defined grid restricts the intervalley scattering to the processes involving energy values up to meV above the conduction band edge. Average deformation potentials of eV/Å for -TA and eV/Å for -LA processes have been obtained.Tyuterev:2010 () These values are close to the values of eV/Å and eV/Å used in Monte-Carlo simulations (column of Table 1).Jacoboni:1983 ()
In contrast with the averaged values for the -TA and -LA processes however, our ab initio averaged values for the intervalley scattering by other phonons can be widely different from the set of deformation potentials of Ref. Jacoboni:1983, . This appears to be especially true for the scattering by the -LO phonon (Table 1), and illustrates the need for parameter free calculations. Moreover, our calculations provide additional values for the scattering involving -TO or -LO phonons, which have been neglected so far in the modeling based on empirical deformation potentials (Table 1). Last but not least, ab initio calculations enable us to discriminate between the deformation potentials of the transverse phonons, which have been considered as degenerate in previous works.Jacoboni:1983 () While this approximation is valid for the -TA intervalley processes, because the ab initio values of -TA and -TA are very similar, some caution must be taken with the scattering involving -TO or -TO phonons, for which we have found substantially different values of the deformation potentials (Table 1).
Iv Results and discussion
To make a valid comparison with the experimental results, the effects of all of the scattering mechanisms described in section II have been taken into account. In Fig. 1, the electronic mobility has been computed with both sets of intervalley deformation potentials presented in Table 1, i.e. those computed ab initio (solid line) and those of Jacoboni and Reggiani, Ref. Jacoboni:1983, (dashed line). The parameters for the other scattering processes were taken from Ref. Fischetti:1991, . In both cases we have used the DFT band structure. The comparison between our calculated electronic mobility and the experimental one (symbols) shows an excellent agreement over the whole range of carrier concentrations.
The mobility calculated using the empirical deformation potential values from Ref. Fischetti:1991, is higher than that obtained from the DFT-calculated ones. In fact the empirical deformation potentials of Jacoboni and Reggiani, Ref. Jacoboni:1983, , were adjusted to obtain the best fit to experimental data. We attribute the difference between the two curves to the increase in the number of scattering channels for the intervalley scattering computed ab initio (see Table 1). In the high concentration region (cm), the electron-phonon interaction is negligible, and the two curves coincide.
Comparing to the theoretical work of Ref. Restrepo:2009, (Fig. 1, crosses), which used from first principles calculations, our results give somewhat better agreement with the experimental data. To our knowledge, the reason does not come from a difference in the treatment of the electron-phonon scattering, as the RTA approximation applied to equation (2) of Ref. Restrepo:2009, amounts to the calculation of the average deformation potentials, as we have done in Table 1. Rather, several reasons can give rise to the observed differences. Firstly, the solution of the BTE was performed in the RTA approximation in the work Restrepo:2009, , i.e. only matrix has been taken into account in eq. 3. We have improved over the RTA approximation by taking into account . Secondly, the elastic scattering by ionized impurities has been computed ab initio in Ref. Restrepo:2009, , whereas we still rely on the model of Jacoboni and ReggianiJacoboni:1983 (), adjusted by the best fit to the experimental data. Lastly, in contrast to work of Ref. Restrepo:2009, , we have included electron-plasmon interaction as in Ref. Fischetti:1991, , which becomes important in the high doping regime (10 cm).
Perfect agreement of our results with the existing experimental dataGeballe:1955 () is shown in Fig. 2 for the Seebeck coefficient . No significant difference in has been found between the values calculated using the ab initio deformation potentials and those obtained with the empirical parameters of Ref. Jacoboni:1983, .
Thermal and electric transport properties of Si have been intensively studied in the past.Jacoboni:1983 (); Fischetti:1991 () However, very few works gave thermoelectric coefficients over a wide range of doping concentrations and temperatures in a systematic way. In order to fill this gap, we have systematically computed electronic mobility , electrical conductivity , Seebeck coefficient and thermal conductivity as a function of the temperature for different doping concentrations, which correspond to different values of the chemical potential . In Fig.3 (a), significant increase of can be observed when the temperature decreases below K, in agreement with the inverse power scaling law , where is a positive number increasing with .Ziman:1960 () The electrical conductivity increases with , since more carriers are generated at higher temperatures for a given , despite the decrease of with temperature. It can be seen that, for a given , has the highest value when the value of is close to the bottom of the lowest conduction band (panel b). The value of , the thermal power, decreases when the temperature increases at a given chemical potential (panel c), as expected from eq. 9. The thermal conductivity increases with the temperature at a given . This trend is however different from that observed for nanowires, where thermal transport is dominated by phonons.zhao2010 ()
We have also studied the effect of doping on the thermoelectric properties at given temperatures. The value of the electronic mobility is almost constant in the low concentration region of cm, and then it rapidly decreases when becomes larger than cm (panel (a) of Fig.4). It can be seen that linearly increases with increasing for weak doping, since is constant in this doping range (panel (b) of Fig.4). At higher carrier concentrations, the increase becomes less significant due to the decrease in the mobility. Furthermore, decreases almost linearly with increasing (panel c of Fig.4). The electronic contribution to the thermal conductivity also sensitively depends on the carrier density in the conduction band (panel (d) of Fig.4). We can see that is very weak in the weak doping domain (cm), as significant thermal transport by electrons can be observed only at high doping levels.
In conclusion, we have combined the Boltzmann transport equation with from first principles calculations of the electronic band structure and of the electron-phonon coupling constants for the intervalley scattering in silicon. This approach is developed to improve the calculation accuracy beyond the relaxation time approximation for the determination of thermoelectric coefficients of semiconductors. The computed electronic mobility and the Seebeck coefficient at room temperature have been compared to experimental data. Good quantitative agreement has been obtained. Furthermore, temperature and doping effects on thermoelectric coefficients have been investigated in a systematic way, and provide predictive data.
Results have been obtained with the (modified) Quantum Espresso package.Giannozzi:2009 (); Baroni:2001 (); Pwscf () We thank Paola Gava for the careful reading of the manuscript and many suggestions for its improvement. Useful discussions with M. Calandra, M. Lazzeri and F. Mauri are acknowledged. The authors acknowledge support from the ANR (project PNANO ACCATTONE), and computer time granted by GENCI (project 2210).
- (1) A. Boukai, Y. Bunimovich, J. Tahir-Kheli, J.-K. Yu, W. Goddard III, and J. Heath, Nature 451, 168 (2008)
- (2) T. Harman, P. Taylor, M. Walsh, and B. LaForge, Science 297, 2229 (2002)
- (3) A. Hochbaum, R. Chen, R. Delgado, W. Liang, E. Garnett, M. Najarian, A. Majumdar, and P. Yang, Nature 451, 163 (2008)
- (4) R. Venkatasubramanian, E. Siivola, T. Colpitts, and B. O’Quinn, Nature 413, 597 (2001)
- (5) D. Wang, L. Tang, M. Long, and Z. Shuai, J. Chem. Phys. 131, 224704 (2009)
- (6) M.W. Oh, D.M. Wee, S.D. Park, B.S. Kim, and H.W. Lee, Phys. Rev. B 77, 165119 (2008)
- (7) B.-L. Huang and M. Kaviany, Phys. Rev. B 77, 125209 (2008)
- (8) D.J. Singh, Phys. Rev. B 76, 085110 (2007)
- (9) G.B. Wilson-Short, D.J. Singh, M. Fornari, and M. Suewattana, Phys. Rev. B 75, 035121 (2007)
- (10) H. Hazama, U. Mizutani, and R. Asahi, Phys. Rev. B 73, 115108 (2006)
- (11) G. Madsen, J. Am. Chem. Soc. 128, 12140 (2006)
- (12) G. Madsen and D. Singh, Comput. Phys. Commun. 175, 67 (2006)
- (13) F. Ishii, M. Onoue, and T. Oguchi, Physica B: Condens. Matter 351, 316 (2004)
- (14) T. Thonhauser, T.J. Scheidemantel, J.O. Sofo, J.V. Badding, and G.D. Mahan, Phys. Rev. B 68, 085201 (2003)
- (15) F. Murphy-Armando and S. Fahy, Phys. Rev. B 78, 035202 (2008)
- (16) O. Restrepo, K. Varga, and S. Pantelides, Appl. Phys. Lett. 94, 212103 (2009)
- (17) T. Dziekan, P. Zahn, V. Meded, and S. Mirbt, Phys. Rev. B 75, 195213 (2007)
- (18) D. Yu, Y. Zhang, and F. Liu, Phys. Rev. B 78, 245204 (2008)
- (19) D. Broido, M. Malorny, G. Birner, N. Mingo, and D. Stewart, Appl. Phys. Lett. 91, 231922 (2007)
- (20) A. Ward, D.A. Broido, D.A. Stewart, and G. Deinzer, Phys. Rev. B 80, 125203 (2009)
- (21) C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645 (1983)
- (22) M.V. Fischetti, Phys. Rev. B 44, 5527 (1991)
- (23) D. L. Rode, phys. stat. sol. (b) 53, 245 (1972)
- (24) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt-Saunders, Philadelphia, 1976) pp. 244–262
- (25) D. L. Rode, “Semiconductors and semimetals,” (Academic Press, New York, 1995) Chap. Low-field electron transport, p. 1
- (26) B. R. Nag, Theory of Electrical Transport in Semiconductors (Pergamon Press, Oxford, 1972) pp. 135–138
- (27) J. Sjakste, N. Vast, and V. Tyuterev, Phys. Rev. Lett. 99, 236405 (2007)
- (28) J. Sjakste, V. Tyuterev, and N. Vast, Phys. Rev. B 74, 235216 (2006)
- (29) J. Sjakste, V. Tyuterev, and N. Vast, Appl. Phys. A 86, 301 (2007)
- (30) J. Sjakste, N. Vast, and V. Tyuterev, J. Lumin. 128, 1004 (2008)
- (31) V. Tyuterev, J. Sjakste, and N. Vast, Phys. Rev. B 81, 245212 (2010)
- (32) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964)
- (33) W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965)
- (34) S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58, 1861 (1987)
- (35) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001)
- (36) P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Phys. Rev. B 43, 7231 (1991)
- (37) in Landolt-Börnstein, Semiconductors: Physics of Group IV Elements and III-V Compounds, Vol. 17a, edited by O. Madelung, H. Weiss, and M. Schulz (Springer-Verlag, Berlin, 1982)
- (38) R. Brent, Algorithms for Minimization without Derivatives (Prentice Hall, 1972)
- (39) G. Gilat and L. J. Raubenheimer, Phys. Rev. 144, 390 (1965)
- (40) The average over states in the final valley only for -processes yields values of the deformation potentials very similar, namely (in eV/Å): =0.60, =1.20, =1.19 and =4.10.
- (41) W. J. Patrick, Sol. Stat. Electr. 9, 203 (1966)
- (42) F. J. Morin and J. P. Maita, Phys. Rev 96, 28 (1954)
- (43) I. Granacher, J. Phys. Chem. Solids 24, 231 (1967)
- (44) T. H. Geballe and G. W. Hull, Phys. Rev. 98, 940 (1955)
- (45) J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids (Oxford University Press, Oxford, 1960) pp. 421–444
- (46) Z. Wang and N. Mingo, Appl. Phys. Lett. 97, 101903 (2010)
- (47) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. De Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. Seitsonen, A. Smogunov, P. Umari, and R. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009)
- (48) S. Baroni, S. de Gironcoli, A. dal Corso, P. Giannozzi, D. Alfé, F. Antoniella, G. Ballabio, C. Bungaro, G. Cipriani, M. Cococcioni, A. Debernardi, G. Deinzer, G. Fratesi, R. Gebauer, M. Lazzeri, K. Mäder, F. Mauri, A. M. Conte, P. Pavone, G. Roma, K. Stokbro, T. Kokalj, and R. Wentzcovitch, http://www.pwscf.org