Temperature Dependence of the Thermal Conductivity of Thin Silicon

pubs.acs.org/NanoLett

Temperature Dependence of the Thermal Conductivity of Thin Silicon Nanowires Davide Donadio,*,† and Giulia Galli†,‡ †

Department of Chemistry and ‡ Department of Physics University of California Davis, One Shields Avenue, Davis, California 95616 ABSTRACT We compute the lattice thermal conductivity (κ) of silicon nanowires as a function of temperature by molecular dynamics simulations. In wires with amorphous surfaces κ may reach values close to that of amorphous silicon and is nearly constant between 200 and 600 K; this behavior is determined by the presence of a majority of nonpropagating vibrational modes. We develop a parameter-free model that accounts for the temperature dependence observed in our simulations and provides a qualitative explanation of recent experiments. KEYWORDS Thermal conductivity, thermoelectric materials, silicon nanowires, molecular dynamics, lattice dynamics

S

dynamics, utilizing the Green function formalism or the solution of the Boltzmann transport equation (BE).15-19 In turn, for core/shell structures, mostly continuum models have been developed.9,20-24 In general, these models relate the decreased lattice thermal conductivity in core/shell nanowires to the diffusive character of the phonon scattering at the crystalline/amorphous interface. With few suitably chosen empirical parameters, they fit well the experimental data for larger diameter SiNW; however, such models are not expected to be applicable for thin wires (d < 10 nm). Nevertheless, even disregarding the atomistic structure of SiNW, the parameter-free model proposed by Martin et al.25 achieves good agreement with experiments for both VLS10 and EE4 SiNW. In a previous work,14 we have shown that the presence of a thin amorphous surface layer may reduce the thermal conductivity of thin SiNW (diameter e 4 nm) by up to 100 times with respect to the bulk value. While ref 14 focuses on the effects of surface disorder on heat transport at room temperature, in this work we compute the thermal conductivity of thin wires as a function of the temperature by equilibrium MD simulations. We then present a parameterfree model based on the BE and a Kubo based theory in the harmonic approximation,26,27 which is in very good agreement with MD simulations in the temperature range between 200 and 600 K. This model allows for the evaluation of quantum effects at low temperature, thus providing a physical interpretation of the unexpected trends of κ(T) at low T, detected experimentally. We use the Tersoff potential28 to model crystalline/ amorphous CS silicon nanowires with diameters (d) of 2 and 3 nm and crystalline SiNW (1.1 < d < 3 nm), grown along the (100) direction. The structural and electronic properties of the crystalline SiNW are reported in refs 29 and 30. CS SiNW models are generated by heating crystalline wires close to their melting point and then quenching them to

emiconducting silicon and germanium nanowires have raised interest as promising building blocks for nanoscale electronic devices,1-3 and, more recently, as efficient thermoelectric materials.4,5 Their size, growth direction, and surface structure may be tuned by selecting different growth protocols and conditions, varying from bottom-up (vapor-liquid-solid growth, VLS) to top-down (electroless etching, EE) approaches,4,6-9 thus modulating their electronic and transport properties. A basic understanding of heat transport in silicon nanowires (SiNW) is highly desirable, either to predict the performance of nanoelectronic elements, or to drive improved fabrication processes of efficient silicon-based thermoelectric devices. To this aim, materials with high power factor and low thermal conductivity (κ) are needed. κ has proved to be very sensitive to the diameter10 and surface structure of SiNW.4 Experiments have shown that size reduction to ∼20 nm may reduce the thermal conductivity by about 1 order of magnitude with respect to the bulk,10 while further reduction may be achieved by modifying the thickness and the roughness of the amorphous oxide that forms spontaneously at the surface.4 The thermal conductivity of these systems, which are effectively crystalline core/amorphous shell (CS) SiNW, displays unusual temperature dependence, with respect to the behavior of κ(T) in crystalline systems. While in semiconducting crystals one expects κ(T) ∼ T3 for T f 0 K and κ(T) ∼ 1/T at higher temperature after a maximum is reached, experiments for SiNW show nearly linear dependence of κ(T) at low T and a broad plateau above ∼150 K that may extend beyond room temperature.9,10 Thermal transport in pristine crystalline wires has been the subject of several theoretical investigations, based either on molecular dynamics (MD) simulations11-14 or on lattice * To whom correspondence should be addressed. E-mail: [email protected]. Received for review: 09/30/2009 Published on Web: 02/17/2010 © 2010 American Chemical Society

847

DOI: 10.1021/nl903268y | Nano Lett. 2010, 10, 847–851

FIGURE 1. Snapshots of the silicon nanowires considered in this study (10 nm supercells). The bonds between atoms belonging to the amorphous shell in core-shell wires are represented by thinner lines.

FIGURE 2. Thermal conductivity of SiNW computed by equilibrium molecular dynamics (MD) for crystalline wires (left), and core-shell wires (right panel). A 1/T trend (lines in left panel) fits well the MD results for crystalline wires, while the thermal conductivity of core-shell wires is constant with respect to the temperature.

room temperature; the thickness (ta) of the amorphous shell is tuned by modulating the annealing time, and a local order parameter is used to identify the atoms belonging to the amorphous region.31 Although we have computed κ for several CS SiNW with different thicknesses of the amorphous shell, here we restrict our attention to explore the temperature dependence of κ of two of them with ta ∼ 0.3 and 0.5 nm, respectively (Figure 1). We consider the temperature range between 200 and 600 K. In analogy with the diffusion coefficients in Brownian motion theory, the thermal conductivity is computed from equilibrium MD simulations by fitting the mean square displacement of the integrated heat flux R by an Einsteintype relation32,33

〈(Rz(t) - Rz(0))2〉 ) 2κ[t + τ(e-t/τ - 1)]

systematic underestimate of the thermal conductivity, due to an insufficient sampling of the low frequency phonon modes.11,35 The results of our equilibrium MD simulations are summarized in Figure 2. In addition to a difference of more than 1 order of magnitude in the thermal conductivity between crystalline and CS SiNWs at room temperature, we observe qualitatively different temperature dependences of κ. As in bulk crystals, in crystalline SiNW κ(T) decreases monotonically as 1/T, while in CS systems the thermal conductivity is nearly constant in the temperature range considered here. This behavior is analogous to bulk amorphous systems above the Debye temperature. (Our results cannot be compared directly to experiments below room temperature, since thermal transport is dominated by quantum effects, that cannot be taken into account in classical MD simulations.) Different temperature trends reflect the different nature of vibrational modes in the two systems, or, in other terms, different dominant scattering mechanisms; κ is proportional to 1/T when normal and Umklapp phononphonon scattering processes are the dominant ones. κ is constant when the main contribution to scattering comes from structural disorder. To understand in detail our MD results, we analyzed several properties of the vibrational modes of SiNWs, such as density of states, participation ratios, group velocities, and mode polarization. In crystalline systems, both acoustic and optical modes have well-defined propagation and polarization directions, and therefore they can be properly treated as phonons. In CS SiNW, only the lowest frequency (ν e 50 cm-1) acoustic branches propagate like phonons and have finite group velocities (v). The participation ratio of such phonons has equal weight on atoms belonging to the crystalline core and the amorphous shell. These modes exhibit the typical character of acoustic modes in cylindrical symmetry;36 namely, one finds two flexure modes with qua-

(1)

where Rz is the component of R along the wire axis (z). In eq 1, τ provides a rough estimate of the average ballistic N εiri. Here εi is propagation time, and R is defined as R ) Σi)1 the energy density of an atom at the lattice site ri. Since the position vector ri is ill-defined in a periodic system, the quantity R is computed as the integral of the heat flux. (In the calculation of the heat flux, we have assigned the threebody terms contribution to the potential energy to the pivoting atoms. However such choice is arbitrary, we have verified that it provides the same numerical results as other possible choices, as also proven in ref 34.) To achieve converged results with this procedure, long MD simulations are required and size convergence must be tested carefully. Every evaluation of κ has been obtained by averaging over 8 simulations at least 16 ns long, on periodically replicated systems with 10 nm supercells (up to ∼3000 atoms) along the main axis of the SiNW. This size has been chosen after testing finite size effects in supercells as long as 100 nm. The use of simulation cells smaller than 10 nm yields a large © 2010 American Chemical Society

848

DOI: 10.1021/nl903268y | Nano Lett. 2010, 10, 847-–851

3, a general quantitative model of thermal transport in SiNW is obtained

dratic dispersion, and two acoustic modes (torsional and longitudinal) with linear dispersion. The modes at higher frequency have v ) 0, even though they are not localized. We find that localized modes appear for ν g 550 cm-1, and they are expected to have marginal to no influence on heat conduction. Our results indicate that an amorphous coating modifies the character of the majority of the vibrational modes in the wire. This effect becomes more significant when the thickness of the amorphous coating is comparable to the diameter of the wire; in thin wires a very thin amorphous shell is indeed sufficient to change drastically their vibrational and heat transport properties. On the other hand, a combined experimental and theoretical study indicates that the thermal conductivity of VLS wires thicker than 32 nm covered with a thin natural amorphous oxide layer may be modeled as in crystalline systems by a Callaway model improved by using bulk phonon dispersion curves.15 The contribution κi to the total thermal conductivity from propagating phonon modes may be expressed in terms of group velocities vi and relaxation times τi, using the single mode relaxation time approximation of the Boltzmann equation (BE)

κi(q) ) Ci(q)vi2(q)τi(q)

κ)

i

πV2 p2ωi2

∑ |〈i|Jz|j〉|2δ(ωi - ωj)

(2)

(3)

j*i

Here 〈i|Jz|j〉 is the z component of the heat flux operator projected on the i and j vibrational eigenstates. The sum of these terms accounts for more than 85% of the thermal conductivity of a-Si at room temperature.39 Using eqs 2 and © 2010 American Chemical Society

q

(4)

j

The sum over i spawns over the propagating phonon modes, while the sum over j accounts for the contribution of diffusive modes. For crystalline wires, eq 4 reduces to the Boltzmann equation in single mode relaxation time, which is sufficient to reproduce the thermal conductivity values computed by MD (see Figure 3a). On the other hand, to quantitatively model the thermal conductivity of CS SiNW, we use the BE only to treat the propagating phonon modes up to 50 cm-1, while the AF model as in eq 3 provides the contribution to κ of diffusive modes. The calculation is not particularly sensitive to the choice of the crossover frequency, as in the range between 50 and 80 cm-1 Di ∼ vi2τi and thermal transport may be properly described by either BE or the AF model. For low frequency modes, it is instead crucial to consider explicitly the contribution of anharmonic effects, through the BE approach. The δ-function in eq 3 is represented by a Lorentzian, which broadening is chosen between 0.05 and 0.2 meV. The variation of κAF as a function of the broadening is represented by the gray-shaded area in Figure 3b. Yet, we notice that it is mostly the diffusivity of the low-frequency modes that is affected by the broadening, so the variation of the results from eq 4 are much less sensitive to the choice of the Lorentzian broadening. The results obtained using eq 4 are in excellent agreement with those of MD simulations in the temperature range 200-600 K, as shown by Figure 3b. The relative contribution of propagating and diffusive modes depends on the diameter of the wire and on the thickness of the amorphous layer. In the CS wire with d ) 2 nm, phonon modes account for half to 60% of the computed value of κ and diffusive, nonpropagating vibrations account for the other half. In the 3 nm CS wire, the contribution of diffusive modes is slightly larger (∼60%). The majority of the BE contribution to κ in these systems is provided by longitudinal acoustic modes. Our model can be extended to T much below room temperature, where the calculation of κ by MD would not be meaningful, due to the quantum nature of vibrations. In eq 4, the temperature dependence of κ comes from the heat capacity and the lifetimes of propagating phonon modes τi ∼ 1/T, while the diffusivity of nonpropagating modes is temperature independent, as anharmonic scattering contributions are neglected. For crystalline SiNWs, for every phonon i we fit the lifetimes as Ri/T, computed at 300 and 500 K, and we use the fitted constants to extrapolate τi(T) at low T. Plugging such lifetimes into eq 2 together with the classical heat capacity Ci ) kB/V, the red curve in Figure 3a is obtained, which simulates the thermal conductivity of a 2 nm SiNW of infinite length. The good agreement with the

Here Ci is the specific heat per unit volume from the Bose-Einstein quantum statistics distribution. vi are computed by finite differences of ω(q) at the Γ point. Lifetimes τi are defined as the relaxation time of the phonon occupation number. This quantity is proportional to the phonon amplitude squared and is computed by projecting the atom displacements in a MD trajectory on the normal modes of the system.37,38 This approach is sufficiently accurate to reproduce MD calculations of κ at room temperature for Carbon nanotubes35 and crystalline SiNW,14 specifically the one with 2 nm considered in this work. Even if nonpropagating modes do not have finite group velocity, they contribute to thermal transport by a diffusion mechanism. In a Kubo-based harmonic theory developed for the study of thermal transport in amorphous materials26,27 (Allen and Feldman theory, AF), a mode diffusivity (Di) is defined to express the contribution κi ) CiDi from such modes

Di )

∑ ∑ Ci(q)vi2(q)τi(q) + ∑ CjDj

849

DOI: 10.1021/nl903268y | Nano Lett. 2010, 10, 847-–851

reduced dimensionality of SiNW, which are effectively quasi one-dimensional systems; for example, linear temperature dependence of ballistic conductance was predicted for carbonnanotubes.40 Amaximumoccurswhenphonon-phonon scattering starts playing a role (T ∼ 50 K) and anharmonic effects dominate in the 1/T high temperature regime. In CS SiNW, the mean free path of the low frequency acoustic phonons at low temperature is limited by interface scattering as well as by the surface roughness, therefore it is temperature independent. As a consequence eq 4 applied to CS SiNW yields an almost linear behavior of κ(T) (κ ∼ T1.2) at low temperature, similar to that observed experimentally9 (Figure 3c). Nevertheless, κ begins to saturate at much lower temperature than in experiments, where κ(T) grows linearly up to 100 K. Such discrepancy is probably due to the very thin diameter of the SiNW considered in this work. It is worth noting that linear κ(T) dependence for T < 100 K, similar to that reported for Silicon nanowires, was predicted for bulk amorphous Silicon by the Kubo based model adopted in this work39 and eventually confirmed by experiments.41 We then argue that the linear temperature dependence of the thermal conductivity of SiNW observed in experiments is due to an interplay between the effect of dimensionality reduction and the presence of nonpropagating vibrational modes similar to those contributing to heat transport in amorphous Silicon, and in hydrogenated a-Si.42 The diffusivity of such nonpropagating vibrational modes, as expressed in eq 3, does not depend on the temperature; therefore the diffusive nature of the majority of the vibrational spectrum is the reason for the broad plateau of κ(T) in thin SiNW with amorphized surfaces. In experiments, the effect becomes more evident the thinner the wires, and the largertheratiobetweenamorphousandcrystallinematerial,9,10 as confirmed by the model proposed by Martin et al.25 Our description cannot be simply extrapolated to higher temperatures as the propagating modes do not exhibit a clear 1/T decay behavior for CS SiNWs, as found for crystalline SiNW. In addition anharmonic effects may become relevant also for diffusive modes, thus limiting the validity of the harmonic approximation in the AF model. In conclusion, we have presented MD simulations of the lattice thermal conductivity of thin crystalline and core-shell (CS) Silicon nanowires (crystalline core and an amorphous surface), as a function of temperature. The analysis of our data using the BE and the AF theory has allowed us to develop a parameter free model that qualitatively accounts for recent experimental data on κ(T) of SiNWs.4,9 We have shown that in CS wires, the presence of a crystalline/ amorphous interface greatly influences the vibrational spectrum, turning vibrational, propagating modes (the majority in crystalline NWs) into diffusive modes, found in a large portion of the vibrational spectrum of a-Si. It is the combined contribution of propagating and diffusive modes that determines the peculiar temperature dependence of the thermal conductivity in SiNW, as observed in recent experiments.

FIGURE 3. Comparison between the equilibrium molecular dynamics results (open circles with error bars) and the model in eq 4 in crystalline (a), and core/shell (b,c) wires with 2 nm diameter. In panel a, the thermal conductivity of crystalline SiNW as obtained from eq 2 is shown for an infinite wire with classical heat capacity (red solid line), finite wires (0.5, 1, and 4 µm) with quantum heat capacity (black dashed-dotted, solid, dashed lines, respectively). In panel b, different models for the thermal conductivity of core-shell wires are compared; namely the solution of the Boltzmann transport equation (BE) in single mode relaxation time approximation (blue open triangles), the Allen and Feldmann (AF) theory (gray shaded area), and the combination of the two models as in eq 4 (red solid circles). The extrapolation to low temperature of the model in eq 4 for a finite wire of 0.5 µm is shown in panel c.

MD data is a consistency check of the extrapolation procedure. We can then use the correct quantum heat capacity and model the thermal conductivity of finite length SiNW (L ) 0.5, 1, and 4 µm) by correcting the phonon lifetimes for the periodic system (τinf) with a boundary scattering term, via the Matthiessen’s rule, so that 1/τ ) 1/τinf + v/L. In the limit of T f 0, this model yields κ ∼ Ta with a ) 1.5 instead of the standard a ) 3 usually observed for bulk crystalline systems. The reason for such temperature trend is in the © 2010 American Chemical Society

850

DOI: 10.1021/nl903268y | Nano Lett. 2010, 10, 847-–851

Acknowledgment. We thank P.B. Allen and I. Savic for useful discussions. Work supported by DOE/BES under grant DOE/BES-DE-FG02-06ER4626.

(19) Markussen, T.; Jauho, A.-P.; Brandbyge, M. Phys Rev B 2009, 79, No. 035415. (20) Mingo, N.; Yang, L. Phys. Rev. B 2003, 68, 245406. (21) Yang, R.; Chen, G.; Dresselhaus, M. S. Nano Lett. 2005, 5, 1111– 1115. (22) Prasher, R. Appl. Phys. Lett. 2006, 89, No. 063121. (23) Prasher, R.; Tong, T.; Majumdar, A. Nano Lett. 2008, 8, 99–103. (24) Murphy, P. G.; Moore, J. E. Phys. Rev. B 2007, 76, 155313. (25) Martin, P.; Aksamija, Z.; Pop, E.; Ravaioli, U. Phys. Rev. Lett. 2009, 102, 125503. (26) Allen, P. B.; Feldman, J. L. Phys. Rev. B 1993, 48, 12581–12588. (27) Feldman, J. L.; Kluge, M. D.; Allen, P. B.; Wooten, F. Phys. Rev. B 1993, 48, 12589–12602. (28) Tersoff, J. Phys. Rev. B 1989, 39, 5566–5568. (29) Vo, T.; Williamson, A. J.; Galli, G. Phys. Rev. B 2006, 74, No. 045116. (30) Vo, T. T. M.; Williamson, A. J.; Lordi, V.; Galli, G. Nano Lett. 2008, 8, 1111–1114. (31) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Phys. Rev. B 1983, 28, 784. (32) Weitz, D. A.; Pine, D. J.; Pursey, P. N.; Tough, R. J. A. Phys. Rev. Lett. 1989, 63, 1747–1750. (33) Zwanzig, R. Annu. Rev. Phys. Chem. 1965, 16, 67–102. (34) Schelling, P. K.; Phillpot, S. R.; Keblinski, P. Phys. Rev. B 2002, 65, 144306. (35) Donadio, D.; Galli, G. Phys. Rev. Lett. 2007, 99, 255502. (36) Thonhauser, T.; Mahan, G. D. Phys. Rev. B 2004, 69, No. 075213. (37) Ladd, A. J. C.; Moran, B.; Hoover, W. G. Phys. Rev. B 1986, 34, 5058–5064. (38) McGaughey, A. J. H.; Kaviany, M. Phys. Rev. B 2004, 69, No. 094303. (39) Feldman, J. L.; Allen, P. B.; Bickham, S. R. Phys. Rev. B 1999, 59, 3551–3559. (40) Mingo, N.; Broido, D. A. Phys. Rev. Lett. 2005, 95, No. 096105. (41) Zink, B. L.; Pietri, R.; Hellman, F. Phys. Rev. Lett. 2006, 96, No. 055902. (42) Liu, X.; Feldman, J. L.; Cahill, D. G.; Crandall, R. S.; Bernstein, N.; Photiadis, D. M.; Mehl, M. J.; Papaconstantopoulos, D. A. Phys. Rev. Lett. 2009, 102, No. 035901.

REFERENCES AND NOTES (1) (2) (3) (4)

(5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18)

Huang, Y.; Duan, X.; Wei, Q.; Lieber, C. M. Science 2001, 291, 630–633. Gudiksen, M. S.; Lauhon, L. J.; Wang, J.; Smith, D. C.; Lieber, C. M. Nature 2002, 415, 617–620. Dong, Y.; Yu, G.; McAlpine, M. C.; Lu, W.; Lieber, C. M. Nano Lett. 2008, 8, 386–391. Hochbaum, A. I.; Chen, R.; Delgado, R. D.; Liang, W.; Garnett, E. C.; Najaran, M.; Majumdar, A.; Yang, P. Nature (London) 2008, 451, 163–167. Boukai, A. I.; Bunimovich, Y.; Tahir-Kheli, J.; Yu, J.-K.; Goddard, W. A., III.; Heath, J. R. Nature (London) 2008, 451, 168–171. Morales, A. M.; Lieber, C. M. Science 1998, 279, 208–211. Cui, Y.; Lauhon, L. J.; Gudiksen, M. S.; Wang, J.; Lieber, C. M. Appl. Phys. Lett. 2001, 78, 2214–2216. Ma, D. D. D.; Lee, C. S.; Au, F. C. K.; Tong, S. Y.; Lee, S. T. Science 2003, 299, 1875–1878. Chen, R.; Hochbaum, A. I.; Murphy, P.; Moore, J.; Yang, P.; Majumdar, A. Phys. Rev. Lett. 2008, 101, 10551. Li, D.; Wu, Y.; Kim, P.; Shi, L.; Yang, P.; Majumdar, A. Appl. Phys. Lett. 2003, 83, 2934–2936. Volz, S. G.; Chen, G. Appl. Phys. Lett. 1999, 75, 2056–2058. Ponomareva, I.; Srivastava, D.; Menon, M. Nano Lett. 2007, 7, 1155–1159. Yang, N.; Zhang, G.; Li, B. Nano Lett. 2008, 8, 276–280. Donadio, D.; Galli, G. Phys. Rev. Lett. 2009, 102, 195901. Mingo, N.; Yang, L.; Li, D.; Majumdar, A. Nano Lett. 2003, 3, 1713–1716. Mingo, N. Phys. Rev. B 2003, 68, 113308. Mingo, N.; Broido, D. A. Phys. Rev. Lett. 2004, 93, 246106. Markussen, T.; Jauho, A. P.; Brandbyge, M. Nano Lett. 2008, 8, 3771–3775.

© 2010 American Chemical Society

851

DOI: 10.1021/nl903268y | Nano Lett. 2010, 10, 847-–851