Letter pubs.acs.org/NanoLett
Prediction of Giant Thermoelectric Efficiency in Crystals with Interlaced Nanostructure Y. S. Puzyrev,*,† X. Shen,† and S. T. Pantelides†,‡,§ †
Department of Physics and Astronomy and ‡Department of Electrical Engineering and Computer Science, Vanderbilt University, Nashville, Tennessee 37235, United States § Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States S Supporting Information *
ABSTRACT: We present a theoretical study of the thermoelectric efficiency of “interlaced crystals”, recently discovered in hexagonal-CuInS2 nanoparticles. Interlaced crystals are I−III−VI2 or II−IV−V2 tetrahedrally bonded compounds. They have a perfect Bravais lattice in which the two cations have an infinite set of possible ordering patterns within the cation sublattice. The material comprises nanoscale interlaced domains and phases with corresponding boundaries. Here we employ density functional theory and large-scale molecular dynamics calculations based on model classical potentials to demonstrate that the phase and domain boundaries are effective phonon scatterers and greatly suppress thermal conductivity. However, the absence of both structural defects and strain in the interlaced material results in a minimal effect on electronic properties. We predict an increase of thermal resistivity of up to 2 orders of magnitude, which makes interlaced crystals an exceptional candidate for thermoelectric applications. KEYWORDS: Thermoelectric, nanoparticles, density functional theory, ternary, phonon scattering
T
scattering,7−9 as in nanocomposite materials.9−15 The size of the scattering centers needs to be somewhere between the phonon and electron mean-free paths to minimize the concomitant decrease in the electrical conductivity. Nevertheless, increasing the density of scattering centers reduces the electron mean-free path, degrading the electrical conductivity. In this Letter, we employ model calculations to show that, in recently discovered interlaced crystals,16 the thermal conductivity can be reduced by up to 2 orders of magnitude, while the electrical conductivity remains essentially intact. Interlaced crystals have either a wurtzite or zinc-blende tetrahedrally bonded structure and two cations sharing the cation sublattice. Theory and atomic-resolution microscopy (Figure 1) demonstrated that there exist many ordering patterns for the cations with essentially equal energy because they all satisfy the valence-octet rule (in an ABX2 compound, every X atom has two A and two B neighbors, ensuring eight electrons for the four X-A and X-B bonds). Domain and phase boundaries cost no energy and induce no strain.16 The underlying Bravais lattice, shown in Figure 1b remains intact. Since electron transport is not sensitive to atomic masses, one expects that the electronic conductivity is not affected by the boundaries, whereas the mass disparity across boundaries should have a major effect on the thermal conductivity.17 Here we report first-principles density func-
hermoelectric materials facilitate conversion of heat into electrical power in energy-saving technologies. One of the applications is low-energy computing, where small-scale thermoelectric components can become important. After decades of research,1−4 many proposed thermoelectric materials have complicated structures, expensive manufacturing technology, and relatively low efficiency.5,6 The maximum efficiency of a thermocouple is expressed as2 ηmax =
Th − Tc Th
1 + ZTm − 1 1 + ZTm + Tc /Th
(1)
where Tc, Th, and Tm are low, high, and average temperatures, and ZT is the dimensionless figure of merit of a thermoelectric material. ZT is written as
ZT =
S2σ T κ
(2)
where S is the Seebeck coefficient, σ is the electron conductivity, and κ is the thermal conductivity. Thus, the critical parameters characterizing a material’s thermoelectric efficiency is the Seebeck coefficient, which is controlled by the material’s electrical conductivity, and the ratio of its electrical and thermal conductivities, which are largely controlled by the scattering of electron and phonons, respectively, by departures from perfect crystalline order. The requirement of lowering the thermal conductivity is typically achieved by controlling the structure of a thermoelectric material at the nanoscale to enhance phonon © XXXX American Chemical Society
Received: August 12, 2015 Revised: November 19, 2015
A
DOI: 10.1021/acs.nanolett.5b03220 Nano Lett. XXXX, XXX, XXX−XXX
Letter
Nano Letters
(DFT) calculations for a DB2 domain boundary in a large supercell (a domain boundary between two orientations of the WZ2 phase, designated as WZ2/WZ2*; two oppositely oriented domain boundaries must be included in a supercell in order to repeat it periodically; a total of 96 atoms are used in a supercell). We used the Heyd−Scuseria−Ernzerhof (HSE) version of hybrid exchange-correlation functional,20 PAW potentials,21 and plane-wave basis as implemented in the VASP code.22 The screening parameter of the HSE functional was adjusted to 0.09 Å so that the calculation reproduces the experimental band gap of chalcopyrite phase of CuInS2 (1.54 eV).23 Spin−orbit interaction was included. The plane-wave cutoff was set at 368.6 eV. The electrostatic potential across a DB2 is shown in Figure 2. It is clear that electrons crossing a domain or phase boundary Figure 1. Structure of interlaced crystals along the [001] direction with ordered cation sublattices. Z-contrast STEM-ADF image of CuInS2 nanocrystal (from ref 16) shows (a) domain boundaries, marked by black lines and phase boundary marked by a red line; (b) perfect Bravais lattice of the interlaced crystal; (c) crystal structure of phases WZ2 and (d) WZ4. Dotted lines show primitive unit cells, and red arrows show cell vectors for repeating Cu−In units.
tional theory (DFT) calculations of the electronic properties of CuInS2 as a prototype interlaced crystal and show that interlacing indeed has a minimal effect on electronic properties and hence the Seebeck coefficient, which is controlled by the electrical conductivity. DFT calculations are not, however, currently possible for thermal transport calculations, and reliable classical potentials are not generally available for ternary or quaternary materials. We, therefore, report simulations of thermal transport in an interlaced model crystal using a generic set of classical potentials and find that the controlling factor is the mass ratio of the two cations. We predict the possibility of a giant enhancement of the thermoelectric efficiency by the presence of interlacing. Materials in the ABX2 family and the more general ACDX4 family, where C and D cations share the B sites in ABX2, retaining the octet rule, are known to have ZT ranging from 0.0001 in CuInS2 to 0.01 in Mn-doped CuInSe2,18 up to 0.5 in Cu2+xCd1−xSnSe4 compounds.19 Interlacing has not been established in these materials so far. Thus, controlled growth of thin films of the large ABX2/ACDX4 family of materials may yield interlacing. Coupled with efficient doping to enhance electrical conductivity, especially using large-atomic-mass dopants that would further lower the thermal conductivity, may lead to giant ZT. We first address the effect of domain and phase boundaries in interlaced CuInS2 on the electrical conductivity. There are two possible effects: scattering by a discontinuity in the crystalline potential and the possible change in the electron effective mass (caused by a possible anisotropy of the effective mass in the phase of a domain boundary or a discontinuity in the effective mass at a phase boundary). We note that the discontinuity in the crystalline potential at either a domain or phase boundaries arises from the fact that an electron traversing a boundary suddenly encounters a sheet of In atoms at sites that were expected to be Cu sites (and/or the other way around). Immediately after that, a new crystalline periodicity is established. In order to demonstrate the resulting discontinuity in the electrostatic potential we carried out density functional theory
Figure 2. DFT calculation of the electrostatic potential across a domain boundary DB2, WZ2/WZ2*, in CuInS2, plotted along the line through Cu and In atoms. The arrow points to the position of the domain boundary.
would not experience any scattering if the Cu and In pseudopotentials were identical. However, we do not have true antisite defects because the octet rule continues to be satisfied, i.e., the perturbation is more like that of an isovalent impurity without any lattice relaxation. Our DFT calculations have in fact found that a domain or phase boundary does not introduce any energy levels in the band gap. Band discontinuities at phase boundaries may provide an estimate of the hindrance to electron transport across the boundary. Such calculations, however, can only be done for an interface between two semi-infinite crystals. The results are not particularly relevant for macroscopic samples with interlaced domains that are of order 1 nm. A more useful criterion is the scattering potential at a phase boundary. We calculated the electrostatic potential across a WZ2/WZ4 phase boundary. The scattering potential, shown in Figure S1, is very weak and short-ranged. We compare it with the scattering of a non-Coulombic impurity (Ge in Si) and find that the phase-boundary scattering potential is much weaker. Thus, Coulomb scattering from dopant impurities, which must be introduced to increase the conductivity of a semiconductor, would have a bigger effect on mobilities than the scattering by the neutral perturbation potential at phase boundaries. Neutral impurities are known to have a negligible effect on carrier mobilities. Indeed, in engineering-level device modeling of mobilities, the effect is not included. We conclude that domain B
DOI: 10.1021/acs.nanolett.5b03220 Nano Lett. XXXX, XXX, XXX−XXX
Letter
Nano Letters
Nonequilibrium molecular dynamics simulations (NEMD) are needed to extract the Kapitza resistance of a boundary. The size of a supercell that contains two domain or phase boundaries, as required to allow periodic boundary conditions, makes the use of DFT-based NEMD impractical. We must use classical potentials, but such potentials for ternary compounds such as CuInS2 are not available. We will, therefore, pursue model calculations as follows. First, we use DFT to calculate the phonon spectra of supercells containing a single phase or a domain boundary. The calculated phonon densities of states are shown in Figure 4.
and phase boundaries in interlaced crystals have a negligible effect on carrier mobilities. To investigate how isotropic is the electron effective mass and how much it differs among different phases, we used DFT calculations in bulk CuInS2 for each of the WZ phases discussed in ref 16. The energy dispersion curve for the bottom of the conduction band, the Γ point, in the WZ2 structure in the x and y directions, as shown in Figure 3, demonstrates that
Figure 3. DFT calculation of the energy dispersion curves at the bottom of the conduction band for CuInS2 along x and y directions. Figure 4. Domain boundaries of (a) DB2 for WZ2/WZ2* and (b) DB4 for WZ4/WZ4* are shown by solid red lines; a2 directions of ordering patterns are shown by dashed red lines; (c) phonon densities of states of the WZ4 order and DB2 boundary (d) phonon densities of states for DB2 and DB4 boundaries.
the electron effective mass is independent of direction, namely, with values of mx* = 0.17me; my* = 0.16me. For all other WZ structures the electron effective mass variation is within 5%. The bottom of the band near the Γ point is lower than any other band by more than 1 eV for both chalcopyrite and wurtzite CuInS2 structures, as reported by Tomic et al.24 Since the majority of electrons that are responsible for conduction have energies within kT of the Fermi level, the effective mass isotropy is adequately described by the symmetry near the bottom of the band near the Γ point. Other potentially interlaced crystals may also have isotropic effective masses near the Γ point. If the band minimum is not at or near the Γ point in the Brillouin zone and effective masses are anisotropic, the interlacing may not be as effective in enhancing the thermal coefficient. The overall conclusion is that, in an interlaced crystal, depending on the details of its band structure, electrons may effectively see a crystal lattice as shown in Figure 1b, i.e., a perfect Bravais lattice. In such case, the electronic conductivity is ultimately determined by the level of doping and the usual limitations imposed by phonon scattering and Coulomb scattering from ionized dopants. We now examine the effect of domain and phase boundaries on the thermal conductivity. It is a well-known fact that grain boundaries in a polycrystalline material cause phonon scattering and decrease the thermal conductivity. The thermal scattering at grain boundaries is described by the Kapitza resistance.17,25,26 Domain and phase boundaries in interlaced crystals also scatter phonons due to mass disparity of the atoms on the cation sublattice. Unlike the case of electron scattering, where the difference in In and Cu pseudopotentials have a negligible effect at domain and phase boundaries, the phonons encounter a major change in the mass arrangement of the atoms at such boundaries since In is significantly heavier that Cu. The phonon scattering strength also varies for different phonon modes and momenta, depending on the boundary’s configuration. Overall, however, the phonon mean free path in an interlaced crystal is limited by the size of the interlaced domains.
The fact that these densities of states are virtually identical suggests that the Cu−S and In−S spring constants are essentially the same (recall that a S atom always has two Cu and two In neighbors). This result stems from the fact that the pseudopotentials of Cu and In are very similar and, by the same token, similar to the pseudopotential of a group-II element like S. In other words, we could use ZnS classical potentials, as long as we use the correct masses for Cu and In. For practical reasons, we shall use GaN classical potentials21 and consider interlaced crystals of the form ABX2, where A is a group-II element and B is a group-IV element. We will show that the heat flow across a domain or phase boundary is controlled by the ratio Cm of the masses of the atoms A and B, i.e., Cm = mA/ mB. In other words, the lowering of the thermal conductivity in an interlaced crystal is proportional to the mass disparity between the cations sharing the cation sublattice. We performed NEMD calculations of the thermal resistance of boundaries using the Muller−Plathe method.27 In this method, the heat flux remains constant and the thermal resistance is obtained from the temperature profile. The heat flux is expressed as m j = 1/2tA ∑ i (υi ,hot 2 − υi ,cold 2) (3) 2 where t is the time of simulation, A is the supercell’s crosssection normal to the heat flow, mi is the mass of the ith atom, and the sum is over exchanges of the velocities υhot and υcold of the atoms from hot and cold regions, respectively. The atomic velocity exchange is performed to maintain a constant flux and obtain a temperature gradient when the system reaches a steady-state regime. In the steady-state regime, the thermal C
DOI: 10.1021/acs.nanolett.5b03220 Nano Lett. XXXX, XXX, XXX−XXX
Letter
Nano Letters conductivity can be obtained by averaging over the heat flux j and the temperature gradient dT/dx,
κ=
⟨j⟩ ⟨dT /dx⟩
shown in Figure 6a. The phase boundary WZ2/WZ4 causes a conductivity decrease of a factor of 4, similar to the domain
(4)
Molecular dynamics simulations were performed on a 24,000 atom supercell with a temperature gradient across an interlaced crystal, which contained two domain or phase boundaries between hot and cold regions, as shown in Figure 5a. The
Figure 6. Temperature profile for (a) domain boundary DB2, shown by upward triangles, and domain boundary DB4, shown by downward triangles. Well-defined temperature drop appears at the position of the domain boundaries. (b) Temperature profile, shown by circles, for the series of domain boundaries DB2, shown as red dashed lines separated by 1 nm on the left to WZ2 domain. Temperature profile, shown by squares, for a random placement of A and B atoms in the cation sublattice in the domain to the right side of the dotted red line. Mass disparity is mA/mB = 2.5.
boundary with the same temperature dependence. The temperature drop itself extends ∼1 nm, corresponding to several lattice constants, unlike the temperature drop for grain boundaries, where it is typically observed across a single lattice constant spacing.17 Since the average size of domains in an interlaced crystal, as seen in Figure 1, is only about 1 nm, i.e., comparable to the extent of the temperature drop across an isolated boundary, we cannot use the above results to obtain the thermal conductivity of the interlaced crystal. In order to accomplish that objective, we constructed a model crystal with parallel DB2 domain boundaries separated by lg ≈ 1 nm as observed in the interlaced nanoparticles.13 Individual domain boundaries are not visible in the temperature profile, shown in Figure 6b. Instead, a large drop across the entire region with densely packed boundaries is observed. From the temperature profile we obtain a factor of 40 decrease of the thermal conductivity at T = 100 K, and a factor of ∼10 at T = 300 K. Since the temperature drop extends over ∼1 nm, i.e., a similar length scale to the size of the domains, and the crystal has both domain and phase boundaries, we have also considered a limiting case where the phase next to WZ2 is constructed by imposing the octet rule with random arrangement of the A and B atoms in the cation sublattice. This structure is a continuous mesh of the WZx, x = 1,...,∞, phases. In this case the conductivity decrease is even larger, shown in Figure 6, reaching a factor of κp/κg ≈ 1/100 at 100 K. If the cation sublattice is totally disordered, however, i.e., the octet rule is no longer satisfied, the electronic conductivity may be seriously impacted. Thus, interlacing is optimal for thermoelectricity. In summary, we used a model ABX2 system based on the experimentally observed CuInS2 nanoparticles and demonstrated that an interlaced crystal exhibits a decrease in thermal conductivity up to 2 orders of magnitude. This effect is due to the mass disparity in the two cations and is universal for materials obeying tetrahedral rule. Combined with the fact that the electrical conductivity is only minimally affected by interlacing, interlaced crystals are potentially exceptional candidates for thermoelectric materials. Additional increase in thermal conductivity may be obtained by quaternary interlaced crystals, i.e., with three different cation masses. We hope that
Figure 5. (a) Schematic representation of the temperature gradient and heat flux in the molecular dynamics simulations. Low temperature is shown in blue, high temperature is shown in red. (b) Crystal structure of a DB2 domain boundary in an ABN2 interlaced crystal. The solid red line marks the position of DB2, and the dashed line shows the direction of ordering along a2. (c) Temperature profile for different mass ratios Cm, mA/mB = 1, 2.5, 10. For the case of mass disparity mA/mB > 1, a temperature discontinuity appears at the position of the domain boundary.
temperature drop across the boundaries depends on the mass ratio Cm and defines the Kapitza conductance G. The Kapitza conductance can be obtained as the ratio of the thermal flux and the temperature drop ΔT, ⟨j⟩ (5) ΔT The Kapitza resistance is defined by R = 1/G. Given the grain conductivity κg and the Kapitza conductance, the thermal conductivity κp of a polycrystalline material with average grain size lg can be written as G=
κp−1 = κg −1 + (lgG)−1
(6)
We considered two different domain boundaries: DB2, the boundary between WZ2 domains, and DB4, the boundary between WZ4 domains, shown in Figure 4a,b. First, we investigated the temperature gradient for an isolated domain boundary, shown in Figure 5b,c. This is accomplished by constructing a system with relatively large distance lg ≈ 10 nm between the domain boundaries. To show the dependence of the temperature profile on mass disparity, we start with a model where the masses are the same, mA = mB (Cm = 1). As expected, there is no temperature drop across the grain boundary (Figure 5d), resulting in R = 1, G = 0. However, for a system with a mass disparity (Cm > 1) the simulations show a temperature drop ΔT across the domain boundary, resulting in a nonzero R. For DB2 and DB4 domain boundaries, the mass disparity on the cation site was set to Cm = 2.5. At T = 100 K, the Kapitza resistances of DB2 and DB4 are similar and the thermal conductance decreases by a factor of 3 for ∼10 nm domains, as D
DOI: 10.1021/acs.nanolett.5b03220 Nano Lett. XXXX, XXX, XXX−XXX
Letter
Nano Letters
(14) Ibáñez, M.; Zamani, R.; Gorsse, S.; Fan, J.; Ortega, S.; Cadavid, D.; Morante, J. R.; Arbiol, J.; Cabot, A. ACS Nano 2013, 7 (3), 2573− 2586. (15) Lesnyak, V.; George, C.; Genovese, A.; Prato, M.; Casu, A.; Ayyappan, S.; Scarpellini, A.; Manna, L. ACS Nano 2014, 8 (8), 8407− 8418. (16) Shen, X.; Hernández-Pagan, E. A.; Zhou, W.; Puzyrev, Y. S.; Idrobo, J.-C.; Macdonald, J. E.; Pennycook, S. J.; Pantelides, S. T. Nat. Commun. 2014, 5, 55431. (17) Maiti, A.; Mahan, G. D.; Pantelides, S. T. Solid State Commun. 1997, 102 (7), 517−521. (18) Zhai, Y.-T.; Chen, S.; Yang, J.-H.; Xiang, H.-J.; Gong, X.-G.; Walsh, A.; Kang, J.; Wei, S.-H. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84 (7), 75213. (19) Liu, F. S.; Zheng, J. X.; Huang, M. J.; He, L. P.; Ao, W. Q.; Pan, F.; Li, J. Q. Sci. Rep. 2014, 4, 5774. (20) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. J. Chem. Phys. 2006, 118 (18), 8207−8215. (21) Kresse, G.; Joubert, D. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59 (3), 1758−1775. (22) Kresse, G.; Furthmüller, J. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54 (16), 11169−11186. (23) Yoshino, K.; Ikari, T.; Shirakata, S.; Miyake, H.; Hiramatsu, K. Appl. Phys. Lett. 2001, 78 (6), 742−744. (24) Tomić, S.; Bernasconi, L.; Searle, B. G.; Harrison, N. M. J. Phys. Chem. C 2014, 118 (26), 14478−14484. (25) Bagri, A.; Kim, S.-P.; Ruoff, R. S.; Shenoy, V. B. Nano Lett. 2011, 11 (9), 3917−3921. (26) Swartz, E. T.; Pohl, R. O. Rev. Mod. Phys. 1989, 61 (3), 605− 668. (27) Müller-Plathe, F. J. Chem. Phys. 1997, 106 (14), 6082−6085.
the present theoretical results will generate interest to fabricate interlaced thin films to test the predictions.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.nanolett.5b03220. Calculations of electrostatic potential across the phase boundary WZ2/WZ4 of interlaced CuInS2 and compare it to the case of Ge substitutional impurity in Si (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Author Contributions
Y.S.P. performed DFT and NEMD calculation. X.S. performed effective mass DFT calculations. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by Department of Energy grant DEFG02-09ER46554 and by the McMinn Endowment at Vanderbilt University. The computations were performed at the National Energy Research Scientific Computing Center, supported by the DOE Office of Science under Contract No. DE-AC02-05CH11231.
■ ■
ABBREVIATIONS NEMD, nonequilibrium molecular dynamics; DFT, density functional theory; DB, domain boundary REFERENCES
(1) Nolas, G. S.; Sharp, J.; Goldsmid, H. Springer: New York, 2001. (2) Rowe, D. M. CRC Press: Boca Raton, FL, 2006; pp 1−1014. (3) Tritt, T. M.; Subramanian, M. A. MRS Bull. 2006, 31 (03), 188− 198. (4) Tritt, T. M. Annu. Rev. Mater. Res. 2011, 41 (1), 433−448. (5) Minnich, A. J.; Dresselhaus, M. S.; Ren, Z. F.; Chen, G. Energy Environ. Sci. 2009, 2 (5), 466−479. (6) Zebarjadi, M.; Esfarjani, K.; Dresselhaus, M. S.; Ren, Z. F.; Chen, G. Energy Environ. Sci. 2012, 5 (1), 5147−5162. (7) Dresselhaus, M. S.; Chen, G.; Tang, M. Y.; Yang, R. G.; Lee, H.; Wang, D. Z.; Ren, Z. F.; Fleurial, J.-P.; Gogna, P. Adv. Mater. 2007, 19 (8), 1043−1053. (8) Hicks, L. D.; Dresselhaus, M. S. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47 (24), 16631−16634. (9) Poudel, B.; Hao, Q.; Ma, Y.; Lan, Y.; Minnich, A.; Yu, B.; Yan, X.; Wang, D.; Muto, A.; Vashaee, D.; Chen, X.; Liu, J.; Dresselhaus, M. S.; Chen, G.; Ren, Z. Science 2008, 320 (5876), 634−638. (10) Joshi, G.; Lee, H.; Lan, Y.; Wang, X.; Zhu, G.; Wang, D.; Gould, R. W.; Cuff, D. C.; Tang, M. Y.; Dresselhaus, M. S.; Chen, G.; Ren, Z. Nano Lett. 2008, 8 (12), 4670−4674. (11) Ma, Y.; Hao, Q.; Poudel, B.; Lan, Y.; Yu, B.; Wang, D.; Chen, G.; Ren, Z. Nano Lett. 2008, 8 (8), 2580−2584. (12) Singh, S.; Ryan, K. M. J. Phys. Chem. Lett. 2015, 6 (16), 3141− 3148. (13) Ibáñez, M.; Zamani, R.; Li, W.; Cadavid, D.; Gorsse, S.; Katcho, N. A.; Shavel, A.; López, A. M.; Morante, J. R.; Arbiol, J.; Cabot, A. Chem. Mater. 2012, 24 (23), 4615−4622. E
DOI: 10.1021/acs.nanolett.5b03220 Nano Lett. XXXX, XXX, XXX−XXX
