Thermal transport size effects in silicon membranes featuring nanopillars as local resonators

Thermal transport size effects in silicon membranes featuring nanopillars as local resonators

Abstract

Silicon membranes patterned by nanometer-scale pillars standing on the surface provide a practical platform for thermal conductivity reduction by resonance hybridization. Using molecular simulations, we investigate the effect of nanopillar size, unit-cell size, and finite-structure size on the net capacity of the local resonators in reducing the thermal conductivity of the base membrane. The results indicate that the thermal conductivity reduction increases as the ratio of the volumetric size of a unit nanopillar to that of the base membrane is increased, and the intensity of this reduction varies with unit-cell size at a rate dependent on the volumetric ratio. Considering sample size, the resonance-induced thermal conductivity drop is shown to increase slightly with the number of unit cells until it would eventually level off.

pacs:

Current affiliation: ]Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena, California 91125, USA

In semiconducting materials, heat is carried mostly by phonons which are quanta of lattice vibrations [Kittel, 1976]. This provides an opportunity to introduce significant changes to the thermal transport properties by direct engineering of the phonon characteristicswhich are shaped primarily by the phonon band structure and the nature of the underlying scattering mechanisms [Chen, 2000; ?]. Recent reviews survey developments in theory, computation, and experiment pertaining to nanoscale thermal transport in a variety of materials and point to the remarkable possibilities for using nanostructuring as a means for phonon engineering [Cahill et al., 2003; ?; ?].

Thermoelectric energy conversion stands to benefit profoundly from the ability to alter the phonon properties by nanostructuring [Chen et al., 2003; ?], as well as by reducing the material dimensionality [Dresselhaus et al., 2007]. Thermoelectric materials, which generate electricity from heat and vice versa, are characterized by a figure of merit defined as , where is the Seebeck coefficient, is the electrical conductivity, is the thermal conductivity (consisting of a lattice component and an electrical component), and is absolute temperature [Rowe, 2005]. One strategy to improve the value of , particularly in semiconductors, is to reduce the lattice thermal conductivity and attempt to do so without negatively affecting and . A promising approach for achieving this goal is to introduce nanoscale local resonators as intrinsic substructures within, or attached to, a host crystalline material [Davis and Hussein, 2014Honarvar and Hussein, 2016]. The emerging system, called nanophononic metamaterial (NPM), exhibits unique properties that are not attainable in conventional nanostructured media such as nanocomposites [Liu et al., 2012; ?] or nanophononic crystals [Yu et al., 2010; ?]. The substructure resonances, which could be numerous for relatively large substructures, may be tuned to couple with all or most of the heat-carrying phonon modes of the underlying host medium. This atomic-scale coupling mechanism is essentially a resonance hybridization between the wavenumber-dependent wave modes of the host medium (phonons) and the wavenumber-independent vibration modes of the local substructure (vibrons). The outcome is significant reductions in the phonon group velocities across roughly the full spectrum of the host medium which, as a result, causes a lowering of the overall lattice thermal conductivity. In Refs. [Davis and Hussein, 2014] and [Honarvar and Hussein, 2016], the host medium is a free-standing silicon membrane and the resonators take the form of silicon nanopillars extruded from one surface. This setup has the advantage that the thermal conductivity reducing elements (the nanopillars) are placed external to the main body of the membrane where in-plane heat transfer takes place, thus providing an environment for minimum electron scattering within the internal domain of the membrane. Using kinetic theory and atomic-scale simulations, a factor of 2 reduction in the thermal conductivity was reported for very small unit cells with membranes less than 5 nm in thickness. A recent molecular dynamics study examined a wider range of geometric dimensions on the same pillared silicon membrane configuration and demonstrated a reduction in the thermal conductivity by roughly a factor of 3 [Wei et al., 2015]. The resonance-hybridization mechanism described above may be used in conjunction with boundary scattering from rough surfaces to lower  [Neogi et al., 2015,Neogi and Donadio, 2015].
For the concept of an NPM to be realized in practice, the nanostructured material’s characteristic length scales, such as the membrane thickness and the nanopillar height and cross-sectional area, need to be at least on the order of a few tens of nanometerswhich is substantially larger than the length scales considered in the previous studies mentioned above. In this Letter, we study the effect of increasing size, in different ways, on the performance of an NPM. While we still consider relatively small models, we seek to establish the fundamental characteristics that govern the relationships between performance and size in order to provide predictive guidelines for much larger systems. We again consider a freestanding silicon membrane with a periodic array of silicon nanopillars standing on the surfacewhich is the same configuration studied in Refs. [Davis and Hussein, 2014], [Honarvar and Hussein, 2016] and [Wei et al., 2015]. Specifically, we consider three size effects: nanopillar size (keeping the membrane thickness and lattice spacing constant), overall unit-cell size, and finite-structure size (which in practice is the sample size and is represented by the number of unit cells existing in a finite structure).
For all our analysis, we consider a unit-cell model consisting of conventional cells (CC) of silicon forming the base (membrane portion) and CC of silicon forming the nanopillar. A silicon CC consists of an eight-atom cube with a side length of nm (Fig. 1a). With this notation, the NPM unit cell has a membrane thickness of and a nanopillar height of . For simplicity, we will limit our attention to nanopillars with a square cross-section, i.e., . For brevity, a labeling convention is adopted whereby the dimensions of a uniform (unpillared) unit cell are represented as (Fig. 1b) and the dimensions of a pillared unit cell are represented as (Fig. 1c).
Equilibrium molecular dynamics (EMD) simulations are used in all studies, and non-equilibrium molecular dynamics (NEMD) simulations are additionally used for the finite-structure investigation. 1 For all simulations, room temperature, K, is assumed and the Stillinger-Weber empirical potential is used to represent the interatomic interactions [Stillinger and Weber, 1985]. Furthermore, we only consider defect-free crystals. In the EMD simulations, the computational domain consists of one unit cell and standard periodic boundary conditions are applied at the in-plane boundaries [Honarvar and Hussein, 2016]. EMD data are used within the Green-Kubo (GK) formalism [Zwanzig, 1965; ?; ?Schelling et al., 2002] to obtain the lattice thermal conductivity, which we will refer to as from now on for brevity. 2 In contrast to EMD, NEMD predicts by direct application of the Fourier’s law of heat conduction, defined as where is the heat current, is the cross section area (equal to or ), and is the temperature gradient [Schelling et al., 2002Tenenbaum et al., 1982Sellan et al., 2010]. A finite number of unit cells is used in the NEMD simulations. Denoting and as the number of NPM unit cells in the - and -directions, respectively, we consider vs. while keeping . Periodic boundary conditions are applied along the -direction and Langevin heat baths are used to apply a temperature gradient across the -direction where the left-end and right-end temperatures are set to K and K, respectively. 3

Figure 1: (a) Silicon conventional cell with the lattice constant , (b) uniform (unpillared) membrane unit cell, and (c) NPM (pillared) unit cell. The - plane corresponds to the (001) plane of a silicon crystal.

We commence by examining the effect of the nanopillar size. Two NPM unit cells are considered with dimensions CC and CC, respectively. The nanopillar-to-base membrane volume fraction is for the unit cell with a small nanopillar and for the unit cell with a large nanopillar. The corresponding thermal conductivity ratio (defined as the NPM thermal conductivity, , divided by the thermal conductivity of a corresponding uniform membrane sized CC, ) is 0.49 and 0.2, respectively. This indicates that the intensity of thermal conductivity reduction increases with , which is expected because a larger nanopillar exhibits more local resonance modes than a smaller one and therefore yields a stronger effect for a given membrane thickness and nanopillar lattice spacing. More local resonances implies more mode hybridizations and consequently a more intense overall group velocity reduction. It is noteworthy that a value corresponds to a factor of 5 reduction in the thermal conductivity, which is the highest to date for an NPM. We also observe that versus drops with a decreasing rate; this is because as the number of resonances increases the hybridization effect eventually reaches saturation.
We now investigate the effect of unit-cell size for each of the configurations considered above. This is done by proportionally expanding all dimensions of the NPM (and corresponding uniform-membrane) unit cells and tracking the thermal conductivity reduction with this uniform increase in size. For the and cases, we respectively consider unit-cell dimensions CC and CC for the following scaling factor values: 1, 2, 3, 4, 5 and 6. This set covers values of membrane thickness ranging from 3.26 nm ( 1) to 19.55 nm ( 6), advancing in increments of 3.26 nm.
The unit-cell size effect results are shown in Fig. 2a where for each of the two cases is plotted against (and ). The thermal conductivity for a corresponding uniform membrane is also plotted to serve as the reference case. 4 The resulting values are given in the inset. We observe, expectedly, to drop significantly with decreasing due to dispersion modification as well as the reduction of phonon mean free paths (MFP) caused by diffuse boundary scattering at the surfaces as a result of their reconstruction at the equilibrium state [Appclbaum et al., 1976; ?]. These effects been studied experimentally in the literature in the context of a thin silicon layer on a substrate [Asheghi et al., 1997; ?], freestanding silicon membranes [Neogi et al., 2015,Cuffe et al., 2015], and also silicon nanowires [Li et al., 2003]. The curves are observed to similarly increase with increasing unit-cell size (due to increasing thickness) until gradual saturation. It is seen, however, that the thermal conductivity of the NPM with large nanopillars grows with unit-cell size at a lower pace than the NPM with small nanopillars. This behavior is desirable as it suggests that a higher starting value leads to a less negative effect of increasing size. These trends will continue until the unit-cell characteristic length scales exceed the full span of the MFP distributionat this point, nonlinear scattering mechanisms will dominate and the likelihood of occurrences of resonance hybridizations will diminish as a result.

Figure 2: (a) NPM and membrane thermal conductivity as a function of overall unit-cell size. NPM-to-membrane thermal conductivity ratio as a function of size is shown in the inset. (b) Correlation between thermal conductivity reduction (left) and group velocity reduction (right) as a function of size for NPMs exhibiting different values. Central-difference slopes of the curves are plotted in the insets.
Figure 3: (a) Dispersion curves and (b) group velocities for and NPMs for two different unit-cell sizes. The membrane and NPM quantities are represented by dark and light colors, respectively. (c) Corresponding density, , and cumulative, , quantities of the ratio of NPM-to-membrane group velocities: is evaluated for incremental THz windows, and is evaluated for a cumulative window. The dashed and solid lines correspond to and unit-cell sizes, respectively. Density, , and cumulative, , thermal conductivity quantities for are plotted in the background: is normalized with respect to its maximum value, and is normalized with respect to . Low-frequency behavior is highlighted.

.

Figure 4: (a) NEMD simulation setup ( in this schematic) and temperature profile across the finite dimension. (b) Thermal conductivity as a function of sample size, shown in normalized form in the inset. The curved dashed lines show exponential fittings. The horizontal lines represent average values for the EMD results (dashed) and converged values for the NEMD results (solid).

The observed and dependencies in the reductions of the thermal conductivity values are explained by comparing with the corresponding phonon group velocity trends. We conduct lattice dynamics (LD) calculations [Gale and Rohl, 2003] and obtain phonon group velocities for the models considered, as well as two more values, up to . For efficient calculations, the dispersion spectrum is computed for only a low frequency range, THz, using the reduced Bloch mode expansion technique (RBME) [Hussein, 2009]. 5 The group velocities, , where each phonon mode is labeled by wavevector and mode number , are obtained by differentiating the frequency dispersion curves over 129 wavevector points by a 3-point finite-difference scheme. An average group velocity quantity is then obtained for each case, defined as where is the number of wavevector points and is the number of modes between frequencies and . In Fig. 2b, we plot versus and the ratio of the NPM-to-membrane average group velocity also versus . The latter quantity is denoted and is evaluated considering all group velocities from to  THz. The central-difference slopes of the two sets of curves are shown in the insets. We see a clear qualitative correlation between the and curves which indicates that the thermal conductivity trends are driven by the manner by which the resonance hybridizations affect the group velocities for the different values of and . Of particular importance is the frequency distribution of the nanopillar vibration modes with respect to the underlying membrane dispersion curves for the various cases. This comparison is demonstrated more explicitly in Fig. 3, where we show the dispersion and raw group velocity curves for the and cases for unit-cell sizes and . For example, we observe that a sharp contrast between the NPM and membrane group velocities extends over a broader frequency range for compared to which explains why rises with . As for the slopes in Fig. 2b, these drop with because as the number of resonances increases the hybridization-induced changes with size eventually reaches saturation. Figure 3c presents the group-velocity density and cumulative values highlighting strong reductions due to the presence of the nanopillars in the low-frequency range of THz, which as demonstrated by the overlaying thermal conductivity curves 6 is the most dominant range contributing to the overall value of . We note that the low-frequency dominance and the monotonic increases of the difference in the cumulative group velocity curves between the small and large unit-cell sizes up to roughly 3 THz justifies our choices of and  THz for the comparative analysis presented in Fig. 2b.
Finally, we use NEMD simulations to investigate the dependency of the thermal conductivity on the size of a finite structure, which is illustrated in Fig. 4. As the number of unit cells increases, phonons with longer wavelengths become available for carrying the heat and thus the thermal conductivities for the uniform and NPM cases gradually increase until they converge to their large sample size values; see, for example, Refs. [Schelling et al., 2002Tenenbaum et al., 1982Sellan et al., 2010] for length-dependency studies on other material systems. We observe that the difference between the two curves increases with , although slightlypossibly this is because more of the long-wave phonons with relatively high group velocities become available for resonance hybridization as the number of unit cells is increased. That difference is also expected to level off for large . For comparison, corresponding EMD results are shown. While both indicate NPM thermal conductivity reduction, the EMD predictions are higher than their counterparts from the NEMD simulations. This discrepancy between the two sets of predictions is not uncommon because the EMD and NEMD methods depend on different factors for their accuracy [Sellan et al., 2010Landry et al., 2008].
In summary, we have established that the extent of thermal conductivity reduction by resonance hybridization increases with the size of the local resonator but decreases with increased overall unit-cell size. However, the rate of this decrease in performance with unit-cell size is inversely proportional to the volumetric ratio of the resonator (nanopillar) to the host medium (base membrane). Thus an NPM with a high value of will sustain most of its performance in large sizes. We have also shown that sample size has a mild effect on the thermal conductivity reduction. The NPM unit cells investigated range in base-membrane thickness from 3.26 to 19.55 nm. One of the pillared silicon membrane configurations considered exhibit a factor of 5 thermal conductivity reduction compared to a corresponding uniform membrane.
This research was supported by the National Science Foundation (USA) CAREER Grant No. 1254931. The authors are most grateful to Prof. Pierre Deymier and Dr. Nick Swinteck for their assistance with the initial set up and usage of the LAMMPS code for the equilibrium MD simulations. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (USA) (Award No. CNS-0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research (USA).

Footnotes

  1. All MD simulations are performed using the LAMMPS software where the heat flux is evaluated based on a specified stress-based formula [Plimpton, 1995].
  2. For the EMD simulations, the systems are initially equilibrated for 1 ns (with a time step fs) at the specified temperature using the ensemble (zero pressure cell size based on constant number of atoms, pressure and temperature). Then, in the ensemble (constant number of atoms, volume and energy), the simulations are run for an additional 6 ns to collect heat fluxes that are recorded every 4 fs. The 6 ns time span is sufficiently large compared to the longest phonon lifetime. With these parameters, the heat current auto-correlation functions converge within the first 500 ps. The reported thermal conductivities are the average of values from six independent simulations with different initial velocities. Furthermore, the thermal conductivities along both the - and -directions are considered, effectively resulting in an averaging over twelve predicted values.
  3. The NEMD simulations are run using fs for 0.8 ns to reach the steady state and then for another 9.6 ns to obtain an average heat flux and temperature profile; the reported thermal conductivities are the average of two independent simulations.
  4. All uniform membrane calculations are based on CC unit cells, for which the thermal conductivities are converged as a function of the simulation cell size [Honarvar and Hussein, 2016].
  5. RBME is a secondary modal-expansion technique for speeding up band-structure calculations. In the present implementation, a 3-point expansion is conducted whereby eigenvectors are selected at the and X points within the irreducible Brillouin zone to form the reduced basis. This procedure generates dispersion curves at more than an order of magnitude higher speed with errors less than 1% compared to full (non-reduced) calculations.
  6. The frequency-dependent thermal conductivity curves are obtained using the Boltzmann transport equation following the single-mode relaxation time approximation as detailed in Ref. [Davis and Hussein, 2014]

References

  1. C. Kittel, Introduction to Solid State Physics (Wiley New York, 1976).
  2. G. Chen, Int. J. Therm. Sci. 39, 471– (2000).
  3. A. A. Balandin, J. Nanosci. Nanotechnol. 5, 1015 (2005).
  4. D. G. Cahill, W. K. Ford, K. E. Goodson, G. D. Mahan, A. Majumdar, H. J. Maris, R. Merlin,  and S. R. Phillpot, Journal of Applied Physics 93, 793 (2003).
  5. D. G. Cahill, P. V. Braun, G. Chen, D. R. Clarke, S. Fan, K. E. Goodson, P. Keblinski, W. P. King, G. D. Mahan, A. Majumdar, H. J. Maris, S. R. Phillpot, E. Pop,  and S. Li, Appl. Phys. Rev. 1, 011305 (2014).
  6. S. Volz, J. Ordonez-Miranda, A. Shchepetov, M. Prunnila, J. Ahopelto, T. Pezeril, G. Vaudel, V. Gusev, P. Ruello, E. M. Weig, M. Schubert, M. Hettich, M. Grossman, T. Dekorsy, F. Alzina, B. Graczykowski, E. Chavez-Angel, J. S. Reparaz, M. R. Wagner, C. M. Sotomayor-Torres, S. Xiong, S. Neogi,  and D. Donadio, Eur. Phys. J. B 89, 15 (2016).
  7. G. Chen, M. Dresselhaus, G. Dresselhaus, J.-P. Fleurial,  and T. Caillat, Int. Mater. Rev. 48, 45 (2003).
  8. C. J. Vineis, A. Shakouri, A. Majumdar,  and M. G. Kanatzidis, Adv. Mater. 22, 3970 (2010).
  9. M. S. Dresselhaus, G. Chen, M. Y. Tang, R. Yang, H. Lee, D. Wang, Z. Ren, J.-P. Fleurial,  and P. Gogna, Advanced Materials 19, 1043 (2007).
  10. D. M. Rowe, Thermoelectrics handbook: Macro to nano (CRC Press, Boca Raton, Florida, USA, 2005).
  11. B. L. Davis and M. I. Hussein, Phys. Rev. Lett. 112, 055505 (2014).
  12. H. Honarvar and M. I. Hussein, Phys. Rev. B 93, 081412(R) (2016).
  13. W. Liu, X. Yan, G. Chen,  and Z. Ren, Nano Energy 1, 42 (2012).
  14. K. Biswas, J. H. He, I. D. Blum, C.-I. Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid,  and M. G. Kanatzidis, Nature 489, 414 (2012).
  15. J. K. Yu, S. Mitrovic, D. Tham, J. Varghese,  and J. R. Heath, Nat. Nanotechnol. 5, 718 (2010).
  16. B. L. Davis and M. I. Hussein, AIP Advances 1, 041701 (2011).
  17. Z. Wei, J. Yang, K. Bi,  and Y. Chen, J. Appl. Phys. 118, 155103 (2015).
  18. S. Neogi, S. Reparaz, L. F. C. Pereira, B. Graczykowski, M. R. Wagner, M. Sledzinska, A. Shchepetov, M. Prunnila, J. Ahopelto, C. M. S. Torres,  and D. Donadio, ACS Nano 9, 3820 (2015).
  19. S. Neogi and D. Donadio, Eur. Phys. J. B 88, 73 (2015).
  20. All MD simulations are performed using the LAMMPS software where the heat flux is evaluated based on a specified stress-based formula [\rev@citealpnumplimpton1995fast].
  21. F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
  22. R. Zwanzig, Annu. Rev. Phys. Chem. 16, 1798 (1965).
  23. A. J. C. Ladd, B. Moran,  and W. G. Hoover, Phys. Rev. B 34, 5058 (1986).
  24. S. G. Volz and G. Chen, Phys. Rev. B 61, 2651 (2000).
  25. P. K. Schelling, S. R. Phillpot,  and P. Keblinski, Phys. Rev. B 65, 144306 (2002).
  26. For the EMD simulations, the systems are initially equilibrated for 1 ns (with a time step fs) at the specified temperature using the ensemble (zero pressure cell size based on constant number of atoms, pressure and temperature). Then, in the ensemble (constant number of atoms, volume and energy), the simulations are run for an additional 6 ns to collect heat fluxes that are recorded every 4 fs. The 6 ns time span is sufficiently large compared to the longest phonon lifetime. With these parameters, the heat current auto-correlation functions converge within the first 500 ps. The reported thermal conductivities are the average of values from six independent simulations with different initial velocities. Furthermore, the thermal conductivities along both the - and -directions are considered, effectively resulting in an averaging over twelve predicted values.
  27. A. Tenenbaum, G. Ciccotti,  and R. Gallico, Phys. Rev. A 25, 2778 (1982).
  28. D. Sellan, E. Landry, J. Turney, A. McGaughey,  and C. Amon, Phys. Rev. B 81, 214305 (2010).
  29. The NEMD simulations are run using fs for 0.8 ns to reach the steady state and then for another 9.6 ns to obtain an average heat flux and temperature profile; the reported thermal conductivities are the average of two independent simulations.
  30. All uniform membrane calculations are based on CC unit cells, for which the thermal conductivities are converged as a function of the simulation cell size [\rev@citealpnumhonarvar2016PRB].
  31. J. A. Appclbaum, G. A. Baraff,  and D. R. Hamann, Phys. Rev. B 14, 588 (1976).
  32. C. J. Gomes, M. Madrid, J. J. Goicochea,  and C. H. Amon, ASME J. Heat Trans. 128, 1114 (2006).
  33. M. Asheghi, Y. K. Leung, S. S. Wong,  and K. E. Goodson, Appl. Phys. Lett. 71, 1798 (1997).
  34. A. M. Marconnet, M. Asheghi,  and K. E. Goodson, J. Heat Trans. 135, 061601 (2013).
  35. J. Cuffe, J. K. Eliason, A. A. Maznev, K. C. Collins, J. A. Johnson, A. Shchepetov, M. Prunnila, J. Ahopelto, C. M. Sotomayor Torres, G. Chen,  and K. A. Nelson, Phys. Rev. B 91, 245423 (2015).
  36. D. Y. Li, Y. Y. Wu, P. Kim, L. Shi, P. D. Yang,  and A. Majumdar, Appl. Phys. Lett. 83, 2934 (2003).
  37. J. D. Gale and A. L. Rohl, Mol. Simulat. 29, 291 (2003).
  38. M. I. Hussein, Proc. R. Soc. A 465, 2825– (2009).
  39. RBME is a secondary modal-expansion technique for speeding up band-structure calculations. In the present implementation, a 3-point expansion is conducted whereby eigenvectors are selected at the and X points within the irreducible Brillouin zone to form the reduced basis. This procedure generates dispersion curves at more than an order of magnitude higher speed with errors less than 1% compared to full (non-reduced) calculations.
  40. The frequency-dependent thermal conductivity curves are obtained using the Boltzmann transport equation following the single-mode relaxation time approximation as detailed in Ref. [\rev@citealpnumPhysRevLett.112.055505].
  41. E. S. Landry, M. I. Hussein,  and A. J. H. McGaughey, Phys. Rev. B 77, 184302 (2008).
  42. S. Plimpton, J. Comput. Phys. 117, 1 (1995).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
278592
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description