Directional Phonon Suppression Function as a Tool for the Identification of Ultralow Thermal Conductivity Materials

Boundary-engineering in nanostructures has the potential to dramatically impact the development of materials for high- efficiency conversion of thermal energy directly into electricity. In particular, nanostructuring of semiconductors can lead to strong suppression of heat transport with little degradation of electrical conductivity. Although this combination of material properties is promising for thermoelectric materials, it remains largely unexplored. In this work, we introduce a novel concept, the directional phonon suppression function, to unravel boundary-dominated heat transport in unprecedented detail. Using a combination of density functional theory and the Boltzmann transport equation, we compute this quantity for nanoporous silicon materials. We first compute the thermal conductivity for the case with aligned circular pores, confirming a significant thermal transport degradation with respect to the bulk. Then, by analyzing the information on the directionality of phonon suppression in this system, we identify a new structure of rectangular pores with the same porosity that enables a four-fold decrease in thermal transport with respect to the circular pores. Our results illustrate the utility of the directional phonon suppression function, enabling new avenues for systematic thermal conductivity minimization and potentially accelerating the engineering of next-generation thermoelectric devices.

Scientific RepoRts | 7:44379 | DOI: 10.1038/srep44379 thermal conductivity, κ, of a material system composed of a circular pores in a square lattice, finding significant heat transport degradation with respect to the bulk. Then, we use the information provided by S(Λ , Ω) to identify a new structure, based on rectangular pores, that exhibits κ as small as 1 Wm −1 k −1 , well below the amorphous silicon limit (1.8 Wm −1 K −1 ) 19 . As our engineering approach can be applied to any combination of material and geometry, it paves the way to high-throughput search of ultra-low thermal conductivity materials.

Bulk
While the methodology developed in this work is applicable to nanostructures with arbitrary materials and shapes, for the sake of clarity we describe our model specifically in relation to porous Si. The first step for computing heat transport in the presence of nanoscale boundaries is calculation of the thermal conductivity of the corresponding bulk material, which we denote κ bulk . With no loss of generality, throughout the text we consider an isotropic bulk MFP distribution of a thermally isotropic material, K bulk (Λ ), which is related to κ bulk via bulk bulk 0 4 where 〈 f(Ω)〉 4π is the angular average (4π) −1 ∫ 4π f(Ω)dΩ and Ω s ( ) is a versor described in terms of the polar angle φ and the azimuthal angle θ, sin sin sin cos cos (2) As , Eq. 1 reduces to the well-known formula bulk 0 . The term K bulk (Λ ) is calculated via density functional theory 6,7 ; for Si at room temperature, we obtain κ bulk ≈ 155 Wm −1 k −1 , in agreement with previous work 20 .

Macroscopic Limit
In porous materials, the volume removal has a degrading effect on heat flow. Furthermore, if the characteristic length L c of the structure is comparable with the MFPs of heat-carrying phonons, size effects take place and the thermal conductivity is further suppressed. In order to compute the "effective" thermal conductivity κ (which is now a scalar) of an array of aligned pores, we identify a square unit-cell with size L containing a single pore, which is chosen to be circular, and apply a difference of temperature Δ T along x, as shown in Fig. 1(a). Periodic boundary conditions are applied at the boundary of the unit cell. We consider the unit cell sizes L = 10 nm and L = 50 nm, keeping the diameter of the pore fixed so that the porosity of the material is φ = 0.25. Assuming the heat flux, J(r), is known, we use Fourier's law, i.e., is an average along the surface of the hot contact and x is normal to this surface. We note that the thermal conductivity is now a scalar. For structures in which phonon-size effects are negligible, heat transport can be modelled by the heat diffusion equation, described by the heat flux J(r) = − κ bulk ∇ T L (r), where T L (r) is the spatially dependent lattice temperature, computed by the continuity equation ∇ · J(r) = 0. Using Eq. 3, we obtain the general expression for the thermal conductivity reduction, For porous materials with low porosity, macroscopic heat reduction is predicted by the Maxwell-Garnett theory which provides the formula κ/κ bulk = (1 − φ)/(1 + φ) ≈ 93 Wm −1 k −1 21 , in excellent agreement with our finite-volume solver of diffusive heat conduction.

Directional Phonon Suppression Function
In aligned porous materials, the characteristic length L c is roughly the pore-pore distance 22 , which in our cases is much shorter than the MFP of most of the dominant phonons in Si 6 , giving rise to significant phonon size effects. In order to take into account these effects, we employ the recently developed MFP-Boltzmann Transport Equation (MFP-BTE) 23 , L where T(r, Λ , Ω) is the effective temperature of phonons with MFP Λ and direction Ω s ( ) and T L (r) is the effective lattice temperature, given by The weights A(Λ ′ ) ensure energy conservation and are expressed as Once Eq. 5 is solved, we include the dependence on Λ and Ω in Eq. 3 to compute the MFP distribution in the nanostructure, i.e., bulk is the thermal flux in the porous material 23 . After substituting Eqs 5-9 into Eq. 8, we obtain the following relationship where S(Λ , Ω) is the "directional phonon suppression function", given by A hot The directional phonon suppression function S(Λ , Ω) is central to our work and describes the MFP dependence and directionality of phonon suppression caused by the boundaries. Once S(Λ , Ω) is known, the thermal conductivity in the nanostructure is given by bulk 0 4 The physical meaning of S(Λ , Ω) can be better understood if we apply an angular average to both sides of Eq. 10, which gives . The term 〈 S(Λ , Ω)〉 4π is also given by which is the conventional MFP-dependent suppression function. This quantity, also called "boundary scattering", is proven to be effective in MFP-reconstruction experiments and for understanding heat transport regime [25][26][27] . We note that the variable Λ is the bulk MFP; any result of this work can be translated in terms of the MFP in the nanostructure, Λ nano , by means of the transformation Λ nano = Λ 〈 S(Λ , Ω)〉 4π .

Diffusive Limit
In the following section, we derive the diffusive and ballistic limits of S(Λ , Ω) and connect the findings with the heat transport regimes of 〈 S(Λ , Ω)〉 4π . For short MFPs, T(r, Λ , Ω) can be expanded up to first-order spherical harmonics as ref.
28 4 4 Combining with Eq. 5, multiplying both sides by ŝ and applying an angular average, this gives 23,28 In Eq. 17, the second term represents the fact that the medium is non-gray. In fact, under the gray-medium approximation, where the lattice temperature is T L (r, Λ ) = 〈 T(r, Λ , Ω)〉 4π , this term vanishes. In this case, phonons are solved independently, leading to a MFP-dependent lattice temperature. For simplicity and with no loss in generality, we now derive the expression of S(Λ → 0, Ω) for our nanoporous system within the gray approximation. Under these assumptions, Eq. 16 becomes the Laplacian equation ∇ 2 T L (r) = 0, i.e., the standard Fourier's law 29 . By assuming a constant heat flux along the hot contact, the gradient of the lattice temperature is simply where we used =n x. The Cartesian representation of S(Λ → 0, Ω), represented by the surface Ω − Λ → Ω Ω = S r s ( ) ( 0, ) ( ) 0, is plotted in Fig. 2(a). We note two lobes, oriented along x and −x, consistently with the fact that both the direction of the applied temperature gradient and the normal of the cold contact are aligned with x. This symmetry can be also seen by the polar representation 0 which, when applied to Eq. 17, becomes ϕ φ ϕ Λ → = S f ( 0, ) 2 ( )sin 2 , as shown in Fig. 2(b). The diffusive limit of 〈 S(Λ , Ω)〉 4π within the gray approximation can be obtained by simply performing an angular average of the polar suppression function, i.e. which is the value obtained by Fourier's law, as shown in Fig. 1(b). However, by using the general expression for non-gray media from Eq. 14, we obtain a lower value for very small MFPs with a peak around 30 nm. This result suggests that the trend obtained by the full MFP-BTE is due to the interaction among phonons with different MFPs. Specifically, the discrepancy in 〈 S(Λ , Ω)〉 4π for low-MFP phonons between the gray and the non-gray models arises from the fact that diffusive phonons in Eq. 5 tend to thermalize to an effective temperature (plotted in Fig. 1(c)) that also depends on ballistic phonons.

Ballistic Limit
We now investigate the ballistic limit of S(Λ , Ω). For large MFPs, boundary scattering becomes predominant and S(Λ , Ω) starts to depend strongly on the geometry of the material. In this regime, Eq. 5 becomes 5 After combining Eqs 11-21, we have A hot which recovers the well-known behaviour S(Λ → ∞ , Ω) ∝ Λ −1 for the ballistic regime 5 . We note that in this case, the directional suppression function is simply the ratio between the MFP distribution in the nanostructure and that in the bulk, i.e. K nano (Λ → ∞ , Ω) = S(Λ → ∞ , Ω)K bulk (Λ ), and so applies to 〈 S(Λ , Ω)〉 4π . As shown in Fig. 2(c), S(Λ , Ω) is pronounced for φ = 0 and φ = π, whereas it rapidly vanishes for other polar angles. This trend can be explained in terms of the view factor, a geometric parameter that quantifies the possibility of having a direct path between the hot and cold contacts 30 . In porous materials with square pore lattices, most of the heat travels through the spaces between the pores, perpendicular to the applied temperature gradient. The relative contribution of such paths is the view factor. Figure 2(d) superimposes S(Λ , Ω) and the material geometry to better elucidate the relationship between the boundary arrangements and phonon suppression. We also note four sub-lobes corresponding to phonons travelling along directions at 45 degrees with respect to the applied temperature gradient, constituting another set of direct paths. For all the other directions, S(Λ , Ω) describes Multiple Phonon-Boundary (MPB) scattering. In our case, heat transport arising from MPB is negligible, with most of the heat carried by phonons travelling through direct paths.

Material Optmization
Using the insights from above, we identify a new nanopore geometry with the same porosity as the circular nanopore case, that has improved properties. This new geometry consists of rectangular pores in a staggered configuration, as shown in Fig. 3(a). The periodicity is larger than the previous case because of the fixed-porosity requirement. From the flux lines plotted in Fig. 3(b), we note that direct paths are absent; thus, phonons scatter multiple times before reaching the cold contact. For low Kns, S(Λ , Ω) is similar to that of the case with circular pores, with the size of the lobes determined by the porosity function f(φ) of the new configuration. Interestingly, for high Kns, S(Λ , Ω) has six preferred directions, as shown in Fig. 3(d). The amplitudes of these peaks are significantly smaller than those in Fig. 2(d). Remarkably, the computed κ are 4 Wm −1 k −1 and 1 Wm −1 k −1 for L = 50 nm and L = 10 nm, respectively, almost five times smaller than their counterparts with circular pores. We note that, for the case with L = 10 nm, κ is smaller than that of amorphous silicon 19 . Finally, as this very low thermal conductivity is obtained with no change in the porosity, the transport of electrons, which travel diffusively, will not be significantly suppressed. As a result, the configuration with staggered rectangular pores has promise as a high-ZT material.

Conclusions
In summary, we have introduced the directional phonon suppression function, a novel concept that captures the essence of size effects in unprecedented detail. Guided by this new quantity, we have identified a porous structure, based on rectangular pores arranged in a staggered configuration, with ultra-low thermal conductivity and relatively low porosity. Our work furthers the general understanding of boundary-dominated heat transport. In addition, as our approach is suitable for any combinations of material and geometry, this work provides practical guidance for the development of novel, high-efficiency thermoelectric materials.

Method
The cumulative thermal conductivity requires the calculation of the three-phonon scattering time. These calculations were performed with a 32 × 32 × 32 grid in reciprocal space and a 4 × 4 × 4 supercell. Force constants calculations were performed with a 5 × 5 × 5 supercell. The isotope disorder scattering is included in the calculation. Phonon-related calculations were carried out using ShengBTE 20 . Density functional theory calculations were performed with Quantum Espresso 31 . Using a projected augmented wave (PAW) pseudopotential 32 , with a plane-wave cut-off of 320 meV, and a 11 × 11 × 11 k-point mesh. Phonon size effects were calculated using an in-house code. The spatial domain was discretized using the finite-volume method solved over an unstructured grid with approximately 6000 elements. The solid angle was discretized into 48 polar and 12 azimuthal angles. The material was assumed to be infinite along the z-direction. The MFP distribution, assumed to be isotropic, was discretized into 30 MFPs.