Phonon Bottleneck Identification in Disordered Nanoporous Materials

Phonon Bottleneck Identification in Disordered Nanoporous Materials


Nanoporous materials are a promising platform for thermoelectrics in that they offer high thermal conductivity tunability while preserving good electrical properties, a crucial requirement for high-efficiency thermal energy conversion. Understanding the impact of the pore arrangement on thermal transport is pivotal to engineering realistic materials, where pore disorder is unavoidable. Although there has been considerable progress in modeling thermal size effects in nanostructures, it has remained a challenge to screen such materials over a large phase space due to the slow simulation time required for accurate results. We use density functional theory in connection with the Boltzmann transport equation, to perform calculations of thermal conductivity in disordered porous materials. By leveraging graph theory and regressive analysis, we identify the set of pores representing the phonon bottleneck and obtain a descriptor for thermal transport, based on the sum of the pore-pore distances between such pores. This approach provides a simple tool to estimate phonon suppression in realistic porous materials for thermoelectric applications and enhance our understanding of heat transport in disordered materials.


The efficient and inexpensive conversion of heat directly into electricity is a long-sought goal with enormous potential in the clean-energy technology landscape Rowe (1995); Tian et al. (2013). The engineering of thermoelectric materials, however, is particularly challenging because of the interrelation of key physical properties constituting the thermoelectric figure of merit ZT, defined as where is the electrical conductivity, is the lattice thermal conductivity, is the Seebeck coefficient, and the temperature. Nanostructuring offers a powerful way to decouple the electrical and thermal transport. In most semiconductors, the numerator of ZT, also referred to as “power factor,” is maximized at relatively high carrier concentrations, so the dominant electron mean free path (MFP) can be as small as a few nanometers Liao et al. (2015). Conversely, phonons may have much larger MFPs, even on the order of microns Esfarjani et al. (2011). Properly engineered nanostructures are therefore able to scatter phonons more effectively than electrons. Porous materials offer a highly tunable platform thanks to their great degree of structural tunability including pore size, shape, and arrangement, as well as the potential for controllable uniform thin films, high temperature resilience and robust contacts. As an example, the thermal conductivities of nanoporous Si have been measured in many studies with the common finding of a strong suppression of thermal transport, leading to a significant improvement in experimentally measured ZT Hopkins et al. (2011); Song and Chen (2004); Lee et al. (2015); Tang et al. (2010); Yu et al. (2010); Verdier et al. (2017); Perez-Taborda et al. (2016); Lee et al. (2017). On the computational level, several models based on the Boltzmann transport equation (BTE) also have shown low thermal conductivities and revealed significant features of phonon-boundary scattering and fundamental thermal transport in nanoporous materials Hao et al. (2009); Lee et al. (2008); Romano and Grossman (2015). Preliminary attempts aiming at tuning thermal conductivity in nanoporous Si have shown that, even within ordered configurations and with pores of the same size, the pattern of the pores can have a large influence on the resulting thermal transport Romano and Grossman (2014). Although aligned configurations offer a robust platform for controllable experiments, pore disorder is unavoidable, especially at smaller length scales Smith et al. (2016). Recent Monte Carlo calculations Wolf et al. (2014a, b) investigated thermal transport in disordered porous materials with circular pores and concluded that the density of pores along the heat flux direction has a significant influence on thermal conductivity. In this paper, we expand on this concept by developing a method that identifies the actual set of pores representing the highest local resistance to phonon transport. To this end, we use the recently developed first-principles BTE solver Romano et al. (2016) to perform thermal transport calculations in random-pore configurations with pores of circular and square shapes. Then, we establish a correlation between the phonon suppression and the pore arrangement within a given configuration, leading to the identification of the pores constituting the phonon bottleneck. Upon introducing a simple descriptor representing the strength of this collection of pores, we find a correlation between such a parameter and the effective thermal conductivity . This work can be potentially used to estimate the degree of phonon suppression in realistic nanoporous samples while avoiding the computational burden of solving the BTE.

Phonon Boltzmann Transport Equation

Our computational approach is based on our recent implementation of the BTE for phonons, which under the relaxation time approximations, reads as  Romano and Grossman (2015)


where is the bulk MFP distribution, is the effective temperature associated to phonons with MFP and direction, , denoted by the solid angle . The term is a bulk material property, and is an angular average. The RHS of Eq.(1) is the effective lattice temperature, a quantity describing the average phonon energy. The term is obtained by using harmonic and anharmonic forces in connection with density functional theory Broido et al. (2007); Esfarjani et al. (2011). The spatial discretization of Eq.(1) is achieved by the finite-volume (FV) method. The simulation domain is discretized by means of an unstructured mesh, generated by GMSH Geuzaine and Remacle (2009). The phonon BTE requires the solid angle discretization to account for different phonon directions. We use the discrete ordinate method (DOM), a technique that solves the BTE for each phonon direction independently and then combines the solutions by an angular integration Abe (1997). As Si is a nongray material, i.e., has a broad MFP distribution, we need to discretize the MFP space, as well. We reach convergence with 30 MFPs (uniformly distributed in log space) and 576 phonon directions. The algorithm is detailed in Ref. Romano and Di Carlo (2011). The overall solution of Eq. 1 requires solving the BTE thousands of times, leading to an increase in the computational time. However, our solver has been conveniently parallelized and each configuration takes only a few minutes with a cluster of 32 nodes.

The walls of the pores are assumed diffusive, a condition that translates into


where and are the solid angle for incoming and outgoing phonons with respect to the contact with normal . Once Eq.(1) is solved, thermal flux is computed via . The effective thermal conductivity is obtained by using Fourier’s law, i.e., , where = 1 K is the applied temperature and K is the distance between the hot and cold contacts (or the size of the unit cell). To focus on phonon size effects, we normalize the thermal conductivity by its diffusive value, i.e., , where is the thermal conductivity computed by the diffusive heat equation and = 156 WmK is the bulk thermal conductivity Li et al. (2014). Our approach has been validated against experiments on porous silicon Song and Chen (2004); Vega-Flick et al. (2016) and, more recently, on silicon labyrinths Park et al. (2017).

Aligned and Random Pores

We first compute thermal transport in aligned configurations, which we will refer to as “aligned circular” (AC) and “aligned square” (AS). The unit cell comprises a single pore and is a square with size L = 10 nm. Heat flux is enforced by applying a difference of temperature = 1 K along the direction. The porosity is fixed at , and periodic boundary conditions are applied throughout. The computed values for are in both cases around 10 W m K , considerably lower than . The magnitude of heat flux, shown in Fig. LABEL:sub@Fig:10a and in Fig. LABEL:sub@Fig:10b for AC and AS, respectively, indicates that phonon travel mostly near the spaces between pores perpendicular to the applied temperature gradient.

Figure 1: Normalized magnitude of thermal flux for the (a) AC, (b) AS, (c) DC and (d) DS cases. The temperature gradient is imposed along the direction. Phonons prefer to travel in the spaces between the pores, as highlighted by the red areas. In all the configurations the pores arrangement is periodic in both and directions. The blue line represents the phonon bottleneck.

For random-pore (or disordered) configurations, the size of the unit cell is chosen to be L = 40 nm, four times as large as that for the aligned cases, in order to generate significant disorder in the pores arrangement. Sixteen nonoverlapping pores are randomly placed while keeping the porosity fixed to = 0.25, thus allowing a direct comparison with the aligned counterparts. We note that the material is still periodic in that pores crossing the border of the unit cell are repeated in the adjacent unit cells. We compute for two hundred arrangements, one hundred for each shape, which we refer to as “disordered circle” (DC) and “disordered square” (DS). The magnitude of thermal flux for two configurations is shown in Fig. LABEL:sub@Fig:10c and Fig. LABEL:sub@Fig:10d, respectively. We note that the formation of high-flux regions is irregular as it depends on the pore configuration. According to Fig. LABEL:sub@Fig:30a and Fig. LABEL:sub@Fig:30b, respectively, the DC and DS cases are found to have average values 15 and 30 lower than that of their aligned counterparts. Intuitively, the combined effect of small bottleneck and vanishing view factor significantly lowers . In the next section, we will analyze in detail the correlation between the pores arrangement and .

Identification of a Descriptor

In previous work Romano and Grossman (2014), we reported that in nanoporous materials is dictated by the view factor and the pore-pore distances. We note that the view factor is a geometrical feature that describes the ability of a ray to travel across the simulation domain without intercepting the pores Howell et al. (2010). In random-pore configurations, the view factor vanishes because of the disordered pores blocking all the direct paths. It is natural, therefore, to speculate whether the average pore-pore distance in the disordered configurations is correlated with . However, after a regression analysis, we conclude that unlike for the ordered case, such a parameter has only a marginal role for the disordered systems. In fact, rigorously speaking, only the interpore spaces perpendicular to heat flux matter. In order to identify the phonon bottleneck we then analyze the pores configuration in terms of graphs.

To this end, we first compute the pores first-neighbor map, as elaborated in the following. A given periodic configuration has a finite set of pores , where is the number of pores. Given two pores and , we define them to be neighbors if, when moving toward , there is no collision with the surrounding pores. The intersection among polygons is computed by the package PyClipper pyc (). After repeating this procedure for all pore pairs, we obtain a first-neighbor map as shown in Fig. LABEL:sub@Fig:20a. We then build the set of edges , where is the number of edges. Each edge connects two neighbor pores, say and , and points toward increasing coordinates, i.e., is connected to only if , where is the circumcenter of the pore . The resulting graph, , is directed in that its edges are unidirectional. We define a path in as a sequence of vertices such that for , where is the length of the path. An elementary circuit is a path where the only repeating vertexes are the first and the last ones, i.e, . In a complete directed graph, the number of distinct elementary circuits, simply referred to as cycles, is


which grows faster than . Although in our case is not complete, the number of cycles can easily reach a few thousand. Here we identify all possible cycles by using Johnson’s algorithm, which has a time bound  Johnson (1975), where is the number of cycles. As the pores are identified uniquely within the unit cell, every pore shares the same label with its periodic counterpart. Consequently, the first and last nodes of a cycle, although having the same identifier, belong to two different unit cells. For our purposes, we select only cycles whose extreme nodes share the same coordinate, as exemplified in Fig. LABEL:sub@Fig:20b. By doing so, we guarantee that the cycles are perpendicular to heat flow and, therefore, are suitable for the identification of a descriptor of thermal conductivity, as explained in the next section.

Figure 2: (a) Example of a first-neighbor map. Each pore in the unit cell is uniquely labeled. The bottleneck is highlighted by the orange line and, in (b), is represented by an elementary circuit, or cycle.

To identify the bottleneck for each configuration we develop the following algorithm:

  1. For each cycle, , we compute the interpore distances of its constituting pores, . Then, we compute the sum of such distances, i.e., .

  2. From the previous point, we have the set . The bottleneck is then .

Figure 3: Distribution of for the (a) DC and (b) DS cases. The straight, horizontal lines represent the aligned counterparts. Distribution of for the (c) DC and (d) DS cases. The vertical lines refer to the bottleneck for the aligned cases.

The phonon bottleneck is the smallest of the sum of pore-pore distances among all the cycles in a configuration. The effectiveness of in describing nanoscale thermal transport in such structures can be estimated by the Spearman correlation rank (), a quantity that indicates how two variables are monotonically correlated to each other  Lehmann and D’abrera (2006). The first step in computing is ranking the values for and and collecting the result via the vectors and , respectively. Then, we compute


where = 100 is the number of simulations for each shape. For both DC and DS cases, we obtain a significant Spearman correlation (higher than 0.63), suggesting that can be used as a good descriptor. We use this knowledge to understand the distributions for the DC and DS cases in relation to the aligned cases. According to simple geometric considerations, the bottlneck for the aligned cases is simply 17.44 nm and = 20 nm. As shown in Fig. LABEL:sub@Fig:30c, for DC, the average is around ; for DS, almost all the configurations have smaller than that of AS [as shown in Fig. LABEL:sub@Fig:30d], due to the square edges. These results reflect the relative trend in between the aligned and disordered cases, corroborating the use of as a valid descriptor for thermal transport. Moreover, we note that most of the bottlenecks have a number of pores (6-7) higher than that of their aligned counterparts (4). This result confirms that smaller , within configurations with the same porosity, can be achieved with anisotropic pore lattices, where the density of pores is higher along the Cartesian direction orthogonal to the applied temperature gradient Wolf et al. (2014a, b). The introduction of a simple descriptor can be used to estimate the ranking of among different samples with disordered pores, supporting experiments on realistic materials Smith et al. (2016); Tang et al. (2010).


In summary, by performing calculations of thermal transport in disordered porous materials we have quantified the effect of the randomness in pore arrangement on the thermal conductivity. Furthermore, we have devised a method to identify the set of special pores composing the phonon bottleneck, potentially empowering experimentalists with a simple tool to assess thermal conductivity in disordered porous materials.


  1. D. M. Rowe, ed., CRC Handbook of Thermoelectrics (CRC Press, Boca Raton, FL, 1995).
  2. Z. Tian, S. Lee,  and G. Chen, J. Heat Transfer 135, 061605 (2013).
  3. B. Liao, B. Qiu, J. Zhou, S. Huberman, K. Esfarjani,  and G. Chen, Phys. Rev. Lett. 114, 115901 (2015).
  4. K. Esfarjani, G. Chen,  and H. T. Stokes, Phys. Rev. B 84, 085204 (2011).
  5. P. E. Hopkins, C. M. Reinke, M. F. Su, R. H. Olsson, E. A. Shaner, Z. C. Leseman, J. R. Serrano, L. M. Phinney,  and I. El-Kady, Nano Lett. 11, 107 (2011).
  6. D. Song and G. Chen, Appl. Phys. Lett. 84, 687 (2004).
  7. J. Lee, J. Lim,  and P. Yang, Nano Lett. 15, 3273 (2015).
  8. J. Tang, H.-T. Wang, D. H. Lee, M. Fardy, Z. Huo, T. P. Russell,  and P. Yang, Nano Lett. 10, 4279 (2010).
  9. J.-K. Yu, S. Mitrovic, D. Tham, J. Varghese,  and J. R. Heath, Nat. Nanotechnol. 5, 718 (2010).
  10. M. Verdier, R. Anufriev, A. Ramiere, K. Termentzidis,  and D. Lacroix, Physical Review B 95, 205438 (2017).
  11. J. A. Perez-Taborda, M. M. Rojo, J. Maiz, N. Neophytou,  and M. Martin-Gonzalez, Scientific reports 6, 32778 (2016).
  12. J. Lee, W. Lee, G. Wehmeyer, S. Dhuey, D. L. Olynick, S. Cabrini, C. Dames, J. J. Urban,  and P. Yang, Nature communications 8, 14054 (2017).
  13. Q. Hao, G. Chen,  and M.-S. Jeng, J. Appl. Phys. 106, 114321 (2009).
  14. J.-H. Lee, G. A. Galli,  and J. C. Grossman, Nano Lett. 8, 3750 (2008).
  15. G. Romano and J. C. Grossman, J. Heat Transf. 137, 071302 (2015).
  16. G. Romano and J. C. Grossman, Appl. Phys. Lett. 105, 033116 (2014).
  17. B. D. Smith, J. J. Patil, N. Ferralis,  and J. C. Grossman, ACS Appl. Mater. Interfaces 8, 8043 (2016).
  18. S. Wolf, N. Neophytou, Z. Stanojevic,  and H. Kosina, J. Electron. Mater. 43, 3870 (2014a).
  19. S. Wolf, N. Neophytou,  and H. Kosina, J. Appl. Phys. 115, 204306 (2014b).
  20. G. Romano, K. Esfarjani, D. A. Strubbe, D. Broido,  and A. M. Kolpak, Phys. Rev. B 93, 035408 (2016).
  21. D. Broido, M. Malorny, G. Birner, N. Mingo,  and D. Stewart, Appl. Phys. Lett. 91, 231922 (2007).
  22. C. Geuzaine and J.-F. Remacle, Int. J. Numer. Meth. Eng. 79, 1309 (2009).
  23. T. Abe, J. Comput. Phys. 131, 241 (1997).
  24. G. Romano and A. Di Carlo, IEEE Trans. Nanotechnol. 10, 1285 (2011).
  25. W. Li, J. Carrete, N. A. Katcho,  and N. Mingo, Comput. Phys. Commun. 185, 1747 (2014).
  26. A. Vega-Flick, R. A. Duncan, J. K. Eliason, J. Cuffe, J. A. Johnson, J.-P. Peraud, L. Zeng, Z. Lu, A. A. Maznev, E. N. Wang, et al.AIP Adv. 6, 121903 (2016).
  27. W. Park, G. Romano, E. C. Ahn, T. Kodama, J. Park, M. T. Barako, J. Sohn J.and Cho, S. J. Kim, A. M. Marconnet, M. Asheghi, A. M. Kolpak, R. Sinclair,  and K. E. Goodson, Sci. Rep. 7, 6233 (2017).
  28. J. R. Howell, M. P. Menguc,  and R. Siegel, Thermal radiation heat transfer (CRC press, Boca Raton, FL, 2010).
  29. PyClipper 1.0.2. .
  30. D. B. Johnson, SIAM J. Comput. 4, 77 (1975).
  31. E. L. Lehmann and H. D’abrera, Nonparametrics: statistical methods based on ranks. (Springer, New York City, 2006).
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 minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

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 description