Super-Suppression of Long-Wavelength Phonons in Constricted Nanoporous Geometries

In a typical semiconductor material, the majority of the heat is carried by long-wavelength, long-mean-free-path phonons. Nanostructuring strategies to reduce thermal conductivity, a promising direction in the field of thermoelectrics, place scattering centers of size and spatial separation comparable to the mean free paths of the dominant phonons to selectively scatter them. The resultant thermal conductivity is in most cases well predicted using Matthiessen’s rule. In general, however, long-wavelength phonons are not as effectively scattered as the rest of the phonon spectrum. In this work, using large-scale molecular-dynamics simulations, non-equilibrium Green’s function simulations, and Monte Carlo simulations, we show that specific nanoporous geometries that create narrow constrictions in the passage of phonons lead to anticorrelated heat currents in the phonon spectrum. This effect results in super-suppression of long-wavelength phonons due to heat trapping and reductions in the thermal conductivity to values well below those predicted by Matthiessen’s rule.


INTRODUCTION
Heat management and the control of heat flow is a cornerstone for a large number of applications and technologies.Heat generating electronic devices require efficient heat sinks of high thermal conductivity materials, while energy harvesting thermoelectric materials require ultra-low thermal conductivities to increase the efficiency of the conversion process [1].For the latter, nanoporous materials have been very successfully used in many occasions to provide drastic reduction in thermal conductivity [2].Even traditionally inefficient thermoelectric materials like Si, can become highly efficient when nanostructured [3,4], or when nanopores are introduced in the lattice [5].The size of the nanostructured features and the distance between them determines the scattering specifics, typically scattering more effectively phonons of similar mean-free-paths (MFPs) [6,7,8,9].For example, short MFPs are effectively scattered by point defects or alloying, medium MFPs by nanopores, nanoprecipitates or superlattices, and long MFPs by large scale grain boundaries in polycrystalline materials.Additionally, although Mie scattering has been shown to play a role in semiconductor materials due to point defects [10], defects with comparable size to the wavelength of phonons are generally required for (strong) Mie scattering to take place.As a result, deliberately scattering large MFP phonons typically requires the engineering of mesoscale defects, such as grain boundaries, because of the long wavelength nature of large MFP phonons.However, such defects often function as recombination centers and are therefore detrimental to electron transport [11].On the other hand, in nanostructures or nano/polycrystalline materials, the long wavelength, long MFP phonons-which are known to carry a large part of the heat due to their large group velocities and weaker scattering-are typically less suppressed compared to other modes [12].In fact, many computational works, by us (which will be discussed below) and others, show that long-wavelength modes are less affected by material defects compared to the rest of the phonon spectrum.This is the case in typical wave-based phonon description formalisms, such as atomistic molecular dynamics (MD), or non-equilibrium Green's functions [13].Even the prominent grain boundary scattering models, for example in particle-based Monte Carlo (MC) methods, or diffusive versus specular scattering models, are wavevector dependent, with the boundary scattering strength or diffuse scattering preference favoring large wavevectors (short wavelengths), rather than the smaller wavevectors (large wavelengths) [14,15].In Refs.[16,17], using Equilibrium Molecular Dynamics (EMD) simulations, we presented a specially designed nanoporous structure of only a small porosity percentage, which allows strong reductions in thermal conductivity.Reductions of even up to 80% larger compared to what one would have expected from Matthiessen's rule, taking into account the pore geometry and distances, can be achieved.We showed that this effect is a result of anticorrelations that appear in the heat current flux.Essentially, heat can be trapped between the pores in such a way that it annihilates itself and effectively leads to drastic phonon MFP reductions.In this work, we revisit this effect.We provide a summary of that observation and the special nanoporous structure design that allows for that.We then advance our understanding of this effect by focusing specifically on the behavior of long-wavelength phonons.We show that in the geometries we describe, we can super-suppress those more resilient long wavelength, long MFP phonons, again largely beyond what Matthiessen's rule suggests.Such large reductions in thermal conductivity with minimal structural intervention to the material, can be largely beneficial to the field of thermoelectrics, for example, in which case the reduced thermal conductivity, in a material which retains high electrical conductivity, would lead to large performance figures of merit.
The paper is presented as follows: In Section II we provide the background information showing the general trend of weak suppression of long wavelength phonons compared to other modes in the presence of nanoscale defects.In Section III we provide our EMD simulation details and describe the structures that show anticorrelated heat current behavior.We further provide explanations for this phenomenon.In Section IV we perform a wavepacket analysis as well as power spectrum analysis of the heat current autocorrelation function of the heat current, and show the signatures of heat confinement and phonon oscillation in the regions between the pores.In Section V we develop a Monte Carlo model to recreate the anticorrelation effects, and compute the suppression function to investigate how the pores affect the phonon spectrum and allow for super-suppression of mean free paths and thermal conductivity.Finally, in Section VI we conclude.

II. WEAK SUPPRESSION OF LONG WAVELENGTH PHONONS
We start by demonstrating the overall resilience of the long wavelength phonon modes to lattice imperfections.For this, we have simulated a series of simple nanoporous silicon geometries (Fig. 1) using the fully wave-based Non-Equilibrium Green's function (NEGF) formalism.The geometries as shown in the sub-figures contain a single or two pores.The channel has fixed boundary conditions surrounding it, i.e. we simulate a nanowire slab, with open boundary conditions at the left/right sides.A complete study and full set of data is presented in Ref. [13], whereas here we focus on the behavior of the long wavelength phonons.
In the first example (Fig. 1a) we increase the diameter of the pore from D = 2 nm to D = 6 nm and subsequently decrease the 'line-of-sight' [18].In the other two examples in Fig. 1b, and 1c, we change the separation of the two pores in the transverse and longitudinal directions, respectively.We then compare the transmission of the different frequency modes between the pristine material and the porous materials as computed from NEGF.As shown in all three cases in Fig. 1a-c, the transmission is in general suppressed for the larger frequency modes, however, the long-wavelength, low-frequency acoustic modes are not affected as strongly (highlighted region to the left of each sub-figure).The larger resilience of these modes persists, despite the large pores and the very narrow line-of-sight that we enforce in the simulation domain [13].This is in contrast to ray-tracing Monte Carlo solutions of the particle-based Boltzmann Transport Equation (BTE) [14], for example, where the narrow line-of-sight reduces the transmission of all phonon modes indiscriminately (this is shown by the solid-flat lines in the Fig. 1 sub-figures).Note that here we are only considering ideal transport conditions, and ignore any effects of phononphonon scattering or mode mixing.Umklapp scattering typically affects high energy phonons more, thus the transmission of high energy phonons is even suppressed compared to that of acoustic long-wavelength phonons to begin with.However, the main point is that the long-wavelength modes are highly resilient to nanostructuring compared to higher frequency modes.In the latter case the neck size, n, is reduced.(b-c) Pores are placed as shown in the schematics.

Geometries and methods:
The structures we model in equilibrium molecular dynamics (EMD), from here onwards, have simplified geometries consisting of cylindrical pores placed in arrays perpendicular to transport as shown schematically in Fig. 2, and extending to the entire depth of the material.We consider a large enough system, of up to 108 nm in length, allowing over 85% of long wavelength modes to transverse the domain [19].The thermal conductivity within EMD is computed via the Green-Kubo method as: where  is the volume of the periodic simulation cell (the volume of Si plus the volume of the pores),  is the temperature,   is Boltzmann's constant,   () is the heat flux in x-, and 〈  ()  ( + )〉 = () is the averaged heat current autocorrelation function (HCACF) at interval .
For our EMD simulations, we use the LAMMPS software [20], and the Stillinger-Weber (SW) potential [21], which is commonly used for heat transport in silicon by us and others [12,22,23,16,17], as it provides a good estimation for the phonon dispersions, despite overestimating the thermal conductivity to values of ~250 W/mK [24].Simulation domain dimensions of up to ~108 nm in length and ~5.4 nm in the width are employed, while the thickness of the domain was ~5 nm in all cases.We typically simulate the same system up to 20 times in order to average the results and reduce the uncertainty in Green-Kubo [18].The simulations were performed at ~300 K and the transport properties are computed along the x-axis, in the [100] crystal orientation.

Anticorrelation effect:
We now introduce the heat current anticorrelation effect that is central to this paper, and the different geometrical features that control it.We consider four porous geometries as depicted in Fig. 3a.The moving average of their HCACFs and the HCACF cumulative integrals from EMD are shown in Figs.3b and 3c, respectively, with the same color as the corresponding geometries in Fig. 3a.
Typically, the HCACF for Si decays exponentially and monotonically to zero [25] (ignoring the expected statistical noise [18,26]), with the decay rate being a measure of the relaxation times of the propagating phonons.The first geometry of Fig. 3a (cyan) with pore radius 1 nm, as we show in Fig. 3b, matches this behavior.It also shows the normal increase and saturation to the cumulative HCACF for the thermal conductivity in Fig. 3c (cyan line).However, the HCACFs of the rest of the three porous geometries with narrow necks show anomalous behavior.They show an anticorrelation trend, where the HCACFs turn negative and then increases again to approach the x-axis from below.This causes a peak in the cumulative HCACF as shown in Fig. 3c, after which a downward trend is observed, reducing the thermal conductivity of the material (see the second geometry in gold color as a basis, in which we increase the size of the pores slightly compared to the first structure).
The dip in the HCACF and the drop in its cumulative integral can be further increased by increasing the pore diameter and reducing the neck size (green structure and data lines).In this case (green structure) a larger dip in the HCACF is achieved, which enhances the reduction in the thermal conductivity in Fig. 3c.In fact, we find that the AC effect is better controlled by the neck to pore radius ratio than either of the two quantities alone [16,17], as the combination of both quantities reflects phonons more effectively.We also find that beyond a certain point, corresponding to ~80% thermal conductivity reductions, reducing the neck does not enhance the AC effect.Possibly a limit is reached at which point the remaining phonons have mean-free-paths smaller than the pore periodicity see Fig. 2), which will thermalize, before anticorrelations form.This suggests that this drastic reduction in thermal conductivity is a result of super-suppression of the long meanfree-path long-wavelength phonons.The fact that the ratio of the neck to pore radius is what determines the AC, is important because it can allow for this effect to be observed for larger neck/pore sizes, which could be more easily realized [17].We can quantify the reduction in  due to AC by considering the distance of the peak of the cumulative HCACF integral from its final converged value.For example, for the geometry in green in Fig. 3c, this value corresponds to a decrease in the accumulated  of approximately 52%.This is an underestimation of the real reduction in thermal conductivity, since in a geometry of similar porosity, but without the AC effects, the cumulative HCAF will continue to further increase beyond the point of the peak in the AC producing geometry, before it eventually saturates.
The anticorrelation effect can also be controlled in the x-axis, i.e. controlling the correlation length of the phonons that contribute to this by increasing the period of the pore array (orange structure compared to the gold) in the transport direction.This shifts the HCACF deep further to the right, at higher correlation times (Fig. 3b), while the peak of the HCACF cumulative integral, Fig. 3c, moves to the right.The overall thermal conductivity is increased compared to the same structure with shorter periodicity, as expected.
We find that there is a linear correlation between the pore periodicity, d, and when the anticorrelation dip minima occur, and this is related to the sound velocity of phonons travelling in between the pore regions.The fact that the dip positioning is related to the sound velocities of the acoustic phonons, is also an indication that the anticorrelation effect targets primarily those modes.
Finally, we note that, many other works in the literature have reported oscillatory or anticorrelation behavior in the HCACF from Green-Kubo simulations, particularly in systems with more than one element.A few examples are metal-organic frameworks (MOFs) [27,28] and organic crystals [29,30,31]; liquids where convective atomic motion (i.e., mass transport) results in negative HCACF minima [32]; layered materials with large atomic mass differences as in HfB2 [33]; in amorphous materials [34,25], in which case the crystalline porous geometries we study here is long-lived (tens of picoseconds).The anticorrelated heat current fluctuations that Haskins et al. [35] reported in carbon nanotubes (CNTs) is the only study that we have found related to ours, as the authors attributed this behaviour to the elastic scattering reflection of the heat flux at the boundary.

Wavepacket analysis and heat trapping:
We have shown that the AC effect can strongly reduce the thermal conductivity, but with limited intervention in the material compared to the typical many scattering centers and large defect surface area.This effect arises from the efficient entrapment of heat between the pores, where phonons undergo reflection before thermalization, and thus strongly reducing their overall mean-free-path.To further illustrate the heat trapping that leads to the AC effect and super-suppression of the MFPs and thermal conductivity, we have modeled the propagation of phonon wavepackets in porous Si and monitored their trajectory.
Wave packets were constructed by preparing the system with the initial atom displacements and velocities arising from a linear superposition of plane waves weighted by a Gaussian distribution around a localized wavevector.In this scheme the displacement (  ) of the  th atom in the  th unit cell along direction μ (x, y or z) due to a wave packet of phonon mode  is given by [17]: where   is the position vector of the  th unit cell,   is the wave vector of the carrier wave and  o and  are the wavepacket's amplitude and momentum uncertainty.The sum is performed over all modes that are in the same acoustic branch as mode γ and that have wave vectors q that are commensurate with the simulation cell.The term   and   are the eigenvector and angular frequency of mode with wave vector q.When preparing a simulation, the initial displacements were computed using  = 0, and the initial velocities were computed from the time derivative of Eq. ( 2).
We set a wavepacket to propagate through a porous Si slab with pore radii of 2 nm and corresponding 1.4 nm neck, a system that shows AC effects.Heatmaps that show the evolution of the kinetic energy of the wavepackets along the material in Fig. 4 for different values, as indicated in each sub-figure.The x-axis shows the propagating direction in this geometry with the position of the pores located at the vertical white lines, separated by 54 nm apart.The y-axis shows the time evolution of the wavepacket kinetic energy.
Fig. 4: Heatmap of the time evolution of the wavepacket kinetic energies in the nanoporous geometry which shows AC effects.The pore radius is 2 nm and the neck size 1.4 nm.The white lines indicate the pore location in the transport direction.Each plot corresponds to a wavepacket of a transverse acoustic mode centered at wavevector, q, as indicated.The distance between the pores is 54 nm.The amplitudes of wavepackets are tuned to keep the temperature as low as 5 K in all of the simulations.Adapted with permission from Ref. [16].Copyrighted by the American Physical Society.
Figure 4 shows that heat, upon encountering the pores, is partially reflected and partially transmitted.Due to the narrow spacing between the pores, phonons (wavepackets here) get trapped and oscillate in the region between the pores.Heat is trapped for all wavevectors.However, note that the wavepackets are only weakly affected by Umklapp scattering, as we consider a very low temperature, 5K, in these simulations.At room temperature, due to their stronger Umklapp scattering, large wavevector/frequency phonons will thermalize before they propagate fully in the region between the pore stacks.Thus, this oscillation is more relevant to long mean-free-path phonons at elevated temperatures, which will bounce a few times, instead of multiple times as we show below.This trapping of low wavevectors that is observed here for the larger pores/smaller necks, leads to drastic reductions in thermal conductivity.

HCACF power spectrum analysis:
To further illustrate indications of heat trapping and super-suppression of long mean-free-path phonons, we now perform an analysis of the spectrum of the HCACF.Various nanoporous geometries (see Fig. 5a) were simulated with EMD as described above.These geometries included pores-which are drilled through the entirety of the slab (left column) -and voids, which are small cavities within the slab (right column).The geometries of these systems varied from uniform placement to staggered, clustered, and randomized placements.Pore size and distance has also been varied between systems (see Table S1 and S2 in the Supporting information file).The colors of the lines plotted in Fig. 5b and 5c correspond to the similarly colored geometries in Fig. 5a.The blue geometry with AC effect (and at a lesser degree the dark orange one), has a large peak with side peaks in the middle region of the depicted power spectrum, originating from the reverberation of phonons between the pore regions.The pore radii (distances) for the geometries in (a) are: 1 nm (20 nm), 1 nm (20 nm), 1.5 nm (51 nm), 1.5 nm (8 nm), 1.5 nm (N/A), 3.5 nm (46 nm), 1.5 nm (8 nm), 1.5 nm (N/A), 1.5 nm (1.6 nm), and 1.5 nm (1.6 nm).The distances listed as N/A are for the porous systems where distances between pores/voids are randomized.
In Fig. 5b, the thermal conductivity calculated via the Green-Kubo formalism has been plotted for each of the nanoporous systems as well as pristine silicon in black (only for the very small correlation times, as it saturates at much higher values compared to the other geometries).A plot showing the full range of the pristine thermal conductivity and the HCACF is shown in the Supporting information in Fig. S1.The first group of points around a porosity of 2% are the first three pore geometries shown in the left column of Fig. 5a.
Despite the low porosity, a significant reduction in thermal conductivity from the pristine to the porous systems is very evident (the black line saturates around 250 W/mK, see Supporting Information Fig. S1b note it is well known that the Stillinger-Weber potential overestimates the thermal conductivity of Si).What is interesting, however, is that the system in this porosity group with the largest, but fewest pores (i.e., with fewer scatterers per volume) reduces the thermal conductivity the most (dark orange bold line).This can be largely attributed to presence of the AC effect arising from that particular backscattering of phonons, as described above.The next geometry group we consider is in the right column of Fig. 5a, with around 5% porosity, consisting of the four void geometries.The same behavior is witnessed for these in Fig. 5b, where thermal conductivity is reduced.The same anticorrelation effect is present, but to a much lesser degree compared to the pore geometries.As voids are not drilled through the entire slab, there are pockets of silicon space around each void where phonons can continue propagating, resulting in fewer backscattering events and less phonon trapping.Now we go back to the first column for the third porosity group of around 10%, which consists of the last three pore geometries in Fig. 5a.Here we see an even more significant reduction in thermal conductivity in Fig. 5b, depending on the size of pores and the distance separating them.In the dark blue geometry, we see the most pronounced anticorrelation effect and the greatest reduction in thermal conductivity.This occurs despite the presence of fewer pores and thus greater distances between the pores, and expectably a larger characteristic scattering length, defined by the distance between the pores along the transport (x-axis) direction.The insert in Fig. 5c shows these thermal conductivity values plotted against porosity for each system.Notice that drastic reductions in thermal conductivity are achieved with very small porositiestypically much larger porosities are needed in common geometries with no AC effects [14].Therefore, pore size and distribution are not the only effects playing a role in thermal conductivity reduction, as evidenced by the different results at equal porosities, but also the distance between pores, both along and perpendicular to the direction of heat transport.
Figure 5c shows the HCACF power spectrum versus the power spectrum frequency here we focus on the left side of the spectrum, for low frequencies for reasons to be evident further below.The power spectrum of the HCACF of phonons that follow the expected Poisson distribution of lifetimes, typically drops with frequency by f --2 , signaling the behavior of Brownian motion, resulting in finite heat flux (see log-log plot in Fig. 6a below for the bulk-black line and the AC devices in dark blue and dark orange).In general, in the presence of pores, in most cases the entire power spectrum reduces in amplitude, and its amplitude becomes more uniform across different power spectrum frequencies, while still a downward trend is typically retained as the frequency increases.
In the case of structures with AC effects (dark blue and dark orange cases), however, additional peaks arise in the power spectrum; these are more pronounced for the dark blue line with the strongest AC in its HCACF.They would correspond to a resonance (e.g. Gaussian or Lorentzian distributions of lifetimes), built on top of the Poisson distribution.
In contrast to Poisson distributions, such resonances lead to zero net heat flux.It would be another indication of phonon reverberation processes between the pores.Essentially, trapped phonons will undergo backward and forward reflections with no contribution to the heat flux.The rise in intensity in those regions of the power spectrum, indicates that a large portion of the heat is trapped into these reverberate oscillations.The majority of those would be long wavelength phonons modes with mean-free-paths that are long compared to the separation between the pore stacks-long enough to undergo (multiple) reflections, and finally thermalize somewhere in between the pore stacks, suppressing their effective meanfree-path.The additional peak for these geometries is also evident in the log-log plot of the power spectrum in Fig. 6a, deviating from the f -2 behavior, indicating a preferable frequency where the system oscillates.
In Fig. 6b, we provide a rough first order estimate of how the power spectrum frequencies of Fig. 5c can be related to corresponding length scales in the system.For this, we proceeded as follows: We have calculated the average speed of the acoustic phonons with the Stillinger-Weber (SW) potential, that is used in this work, which turned out to be υac∼6033m/s [16].By multiplying the τ -1 variable in the x-axis of the spectra in Fig. 5c with υac -1 we can obtain an estimate of an inverse length l -1 , related to path that acoustic phonons will travel before thermalizing.Figure 6b is a repetition of Fig. 5c, but with the x-axis changed to this length.The different colored lines correspond to the same lines in Fig. 5c and the geometries of Fig. 5a.For the specific region of the peak in the power spectrum this will constitute the characteristic reverberation length scale (as we name the x-axis).We show the spectra for the large lengths, i.e. the ones larger than l~80 nm (from the right side of the x-axis towards the origin at the left side).The geometry which shows the supersuppression of thermal conductivity due to the AC effect (thick blue), has its peak centered around ~225 nm.This indicates a characteristic length scale behavior of phonons in the system.The distance between the pores in the propagating direction is 54 nm, thus this could correspond to four periods of phonon travel (e.g. two back and forth reflections between consecutive pore stacks; or one back/forth reflection in two pore stack regions).At this point we don't offer any further analysis on this, the side peaks (whether they correspond to more/less number of oscillations), or even separate the effect of pores on longitudinal or transverse phonons with different sound velocities, since this this provides only a rough estimate.

Fig. 6: (a)
The spectrum of the HCACF in log-log scale for bulk Si (black line) and the AC geometries of Fig. 5a, following the same coloring scheme.(b) The spectrum of the HCACF, versus a first order estimate of a length scale in the system, indicating the reverberation of phonons.The x-axis shows the inverse of the length scale, with the far right region corresponding to ~80 nm.The system with the largest AC effect (thick blue line), peaks around length ~ 225 nm.

V. PHONON SUPRESSION FUNCTION WITHIN A MONTE CARLO MODEL
To elucidate how AC arises from heat trapping between pores, we have developed a statistical model for a grey population of phonons at equilibrium.The model is described in detail in Ref. [17], but is briefly outlined again here.The instantaneous heat flux in a volume  of crystal at thermal equilibrium with no net temperature gradient arises from the stochastic occupancy of the phonon modes in the volume, so that: where   =   ℏ (  ⋅  ̂), is the energy current along  ̂ due to a single phonon in the i-th phonon mode.Here ℏ is the reduced Plank constant, and   and   are the phonon mode's angular frequency and group velocity.The instantaneous occupancy,   (), in Eq. 3 is the number of phonons exciting the i-th mode at time , and the sum is performed over all vibrational degrees of freedom (phonon modes) in the volume.
The key insight of our model is that, while the instantaneous heat flux depends on the superposition of the occupancy   () in all modes, one can decompose the total autocorrelation function into just the sum of correlation functions of individual phonon modes' excursion from their mean Bose-Einstein occupancy Δ  () =   () − 〈  〉 if one assumes there to be no correlation between the initial and scattered phonon mode(s) in each intrinsic scattering event.This approximation is equivalent to the relaxation time approximation (RTA), and results in an expression for the HCACF as: where 〈  〉 is the -th mode's average phonon lifetime, 〈Δ 2 〉 is the variance of the occupancy, and (, ) is the autocorrelation function of the boxcar function Π 0,θ (), with width , so that: The variance of the occupancy is computed by the summing over fluctuations weighted by the probability ) of finding a mode in state of occupancy , and is found to be proportional to the derivative of the average occupancy with respect to temperature The average of the excitation autocorrelation functions, 〈(, )〉, is performed over the Poisson distribution of   = 1 〈〉  − 〈〉 ⁄ of random waiting times for a stochastic process: Putting Eq.s ( 6) and ( 7) into (4) leads to an expression for the HCACF as: where   is the -th vibrational mode's contribution to the volumetric specific heat.When Eq. ( 8) is used in the Green-Kubo formula in Eq. ( 1), it yields the expression for the thermal conductivity derived from the Boltzmann transport equation using the RTA.
The analysis above relies on decomposing the heat current into independent fluctuation events which are self-correlated, but not cross-correlated with one another.The addition of pores to a crystal introduces extrinsic backscattering of phonons at pore surfaces.
In this process, when a phonon is backscattered from mode  →  ′ , the heat current   ′ from the backscattered phonon is correlated with the heat current   from the incident mode because, even in the case of diffuse scattering where the direction of the scattered phonon is independent of the incident phonon, at least the sign of the heat flux normal to the scattering surface is reversed by backscattering making   and   ′ anticorrelated.In this scenario, the independent self-correlated events that make up the total HCACF are still the period between intrinsic (uncorrelated) scattering events, but the heat current during an event now consists of a sequence of discrete values as the phonon is backscattered between different phonon modes.The HCACF from an individual trajectory with multiple backscattering events sequence,   (), is a linear piecewise function that is straight forward to compute numerically.
We have developed a simple ray-tracing scheme to simulate the trajectories of phonon undergoing defuse elastic backscattering and compute their HCACF.The code performs Monte Carlo sampling over trajectories and lifetimes to compute the average HCACF 〈  ()〉 from an isotropic grey phonon gas.The details of this model are described in a prior work [17]; but here we apply the model to compute the phonon suppression function to determine how pores affect the phonon spectrum's contribution to the total thermal conductivity.We used the model to compute the HCACF for a phonon gas with mean free path, , over a broad sweep of the Knudsen number from the diffusive to ballistic regimes.The Knudsen number is defined as the average of the MFPs of grey phonons over the geometrical distance between the pores,   = /.This was done for both the staggered and aligned pore geometries in the inset of  To compare the MC model against the MD simulations we used AlmaBTE [36] to compute the phonon lifetimes for Si using the 2 nd and 3 rd order interatomic force constants for Si modeled with the Stillinger-Weber empirical potential.The distribution of phonon properties computed on a 30 × 30 × 30 -point Brillouin zone mesh where spherically symmetrized and used to compute the HCACF for a "multi-grey" model of monolithic and nanoporous SW Si using: the results of which are shown in Fig. 7a.For fair comparison to the MD simulations where mode occupancy is classical, the classical values of   were used in Eq. 9.This model has no free tuning parameters, and although the intrinsic phonon lifetimes predicted by AlmaBTE are longer than those observed in MD (which include four and more phonon interactions), its overall match of the multi-grey MC model to the MD HCACF is qualitatively very good.The MC computed HCACF shows clear anticorrelation for aligned pores that is absent for staggered pores, as observed in the MD as well.
Having established the validity of the ray tracing model, we used it to compute the thermal conductivity suppression function (  ) for the phonon gas in these two pore layouts with: where  is the volume fraction of the pores.The suppression function indicates how much the pore geometry suppresses the thermal conductivity.The suppression function for the staggered and aligned pores are shown in Figs.7b and 7c, respectively.The suppression in the diffusive regime (low   ) originates entirely from the reduced number of vibrational degrees of freedom as a result of the porosity and is the same for both geometries as they have the same .However, during the crossover to the ballistic regime the aligned pores exhibit a steeper suppression of the thermal conductivity (lower blue line in the inset of Fig. 7b).To better understand this transition, we fitted the suppression function with a Matthiessen rule (MR) model as: where   = (1 − ) is the suppression coming from the missing material at the pores (with  being the volume fraction of porosity), and  is the effective extrinsic phonon meanfree-path for scattering at pores, with  being a dimensionless geometric factor.For the case of the staggered pores, the suppression is fit well by Eq. 11 with geometric scattering factor  = 0.62 (the black line in Fig. 7b).This factor makes sense, as the pores are staggered and so the typical distance that a phonon must travel to impinge on a pore is on the ∼ 1 2 .
In contrast, the suppression function from the aligned pores is not fit well by Eq. 11 using a single value for  (Fig. 7c).The left-hand portion of the suppression plot (the diffusive regime) can be matched by Eq. 11 using a value of  close to 1 (dashed black line), which makes sense as the pores are aligned.However, to match the right-hand (ballistic) portion of the suppression plot to MR requires using a much smaller geometric factor, indicating that in the ballistic regime the effective mean-free-path for scattering phonons at pores is significantly smaller than the distance between the pores.To account for this crossover in behavior, instead of using a fixed value of  in Eq. 11 as would be expected by Matthiessen's rule, we have used a scale-dependent geometric factor: where   , and   are the dimensionless characteristic scattering lengths in the diffuse and ballistic limits, and   is the Knudsen number at the crossover between these regimes.A good fit to the suppression function from aligned pores was obtained with fitting parameters   = 0.92,   = 0.17 and   = 0.19 which is shown by the solid black line in Fig. 7c.This deviation from Matthiessen's rule and the dramatic reduction of the geometric scattering factor from   to   between the diffusive to the ballistic regimes represents a five-fold super-suppression of the phonon mean free path, and a more than three-fold reduction compared to the staggered case!
In the staggered geometry where the pores are better dispersed, the average line of sight distance between pores along the transport direction is shorter and thus in the diffuse regime staggered pores provide more effective thermal resistance.However, in the ballistic regime we observe strong anticorrelation in the HCACF of the aligned pores, which is not present for the staggered pores.In these MC simulations the backscattering behavior when a phonon collides with the pores is identical in both the aligned and staggered cases, thus the anticorrelation must arise not from individual pores, but from the collective arrangement of pores.We propose that this is due to a shadowing effect from the aligned pores.
The anticorrelation in the heat current along  will be strongest from phonons that collide and are backscattered by surfaces perpendicular to  as the anticorrelation arises only in the component of velocity normal to the scattering surface.In both systems, the MC and MD Green-Kubo simulations are performed at equilibrium so while we call what follows a shadowing effect, the phonon radiance hitting the pore surface is the same at all points around the pore.However, the fetch that the phonons travel across before impinging on the pore is not the same in all directions.In the aligned pores, shadowing by neighboring pores means that phonons striking the pores in the ligature sections between pores can only have traveled a short distance before hitting the pore boundary.The contribution that an individual fluctuation event (phonon) makes to the HCACF, is proportional to the square of its distance of travel.Thus, the shadowing does not change the intensity of phonons striking the pores laterally, but diminishes their contribution to the HCACF.This means that the HCACF is composed mainly from long fetch phonons, which also strike the pores in the direction that gives maximum anticorrelation from backscattering.
To understand which phonons are suppressed by the pores, we plot in Fig. 7d the cumulative conductivity over the spectrum of intrinsic phonon mean-free-paths (MFPs) in Si, as well as that of the two geometries examined.While contributions from phonons with all MFPs' to conductivity is reduced, the cumulative distribution becomes particularly flat for the long MFP modes in the porous geometries, and more prominent in the AC geometry (blue line).To understand the correlation between phonon wavelength and backscattering we have plotted in Fig. 7e and 7f a scatter plot of wave number of phonon modes in Si vs their suppression function from the staggered and aligned pores.The suppression function is the factor by which a modes contribution to conductivity is suppressed.We see that the lowest suppression factors are reached in the case of aligned pores in the low and mid  regions highlighted with the blue oval.The same region in the plot for the staggered pores is empty.

Fig. 1 :
Fig. 1: Normalized Si nanochannel phonon transmission functions for different phonon energies ℏω, extracted from NEGF, indicating the resilience of the long wavelength acoustic phonon transmission to nanostructured defects (left highlighted region).The top row shows schematics of the geometries studied.The bottom row shows the transmission functions normalized to the transmission of the pristine channel.The channel has a length of 100 nm and width 10 nm.The flat-solid lines indicate the normalized transmission of the ray-tracing Monte Carlo solution of the BTE.(a) A single pore is placed in the domain with two diameter cases: D = 2 nm (blue) and D = 6 nm (red).In the latter case the neck size, n, is reduced.(b-c) Pores are placed as shown in the schematics.

Fig. 2 :
Fig.2: Schematic of the considered structure to allow for heat current anticorrelations and ultra-low thermal conductivities.Pores are designed with specific separation periods in the transport direction, diameters and neck sizes.The assumed underlying material is Si.

Fig. 3 :
Fig. 3: The anticorrelation (AC) effect.(a) The four porous Si geometries considered: one that does not show AC (blue), and three with AC. b) The HCAF for the four cases.c) The cumulative HCAF.The data coloring in (b) and (c) correspond to the structure coloring in (a).The pore radii (distances) in (a) are: 1 nm (27 nm), 1.5 nm (27 nm), 2 nm (27 nm) and 1.5 nm (54 nm).
McGaughey and Kaviany attribute this behavior to distinct local environments.In all these examples, the mechanisms responsible for anticorrelations are distinct from the diffusive backscattering of the heat flux in the crystalline close-packed nanoporous geometries presented in this work.For example, in our simulations the material is crystalline, consists of only a single element, and has fixed center of mass, so convective heat flux is negligible and only the virial term contributes to the thermal conductivity.Moreover, in these examples, the HCACF oscillations have sub-picosecond period, whereas the HCACF dip in

Fig. 5 :
Fig. 5: Analysis of the spectrum of the HCACF.(a) The geometries simulated with EMD.(b) The HCACF integral, or cumulative thermal conductivity.The black line indicates the corresponding bulk Si HCACF integral.In bold we show the geometries that have strong AC effect.(c) The low frequency end of the power spectrum of the HCACF.The power spectrum is plotted with a Gaussian smoothing of 100 to remove much of the accompanying noise.The data coloring in (b) and (c) correspond to the structure coloring in (a).Inset of (c): the thermal conductivity of the systems versus porosity.The blue geometry with AC effect (and at a lesser degree the dark orange one), has a large peak with side peaks in the middle region of the depicted power spectrum, originating from the reverberation of phonons between the pore regions.The pore radii (distances) for the geometries in (a) are: 1 nm (20 nm), 1 nm (20 nm), 1.5 nm (51 nm), 1.5 nm (8 nm), 1.5 nm (N/A), 3.5 nm (46 nm), 1.5 nm (8 nm), 1.5 nm (N/A), 1.5 nm (1.6 nm), and 1.5 nm (1.6 nm).The distances listed as N/A are for the porous systems where distances between pores/voids are randomized.
Fig. 7a, with a Gaussian process regression used to construct an interpolation function  ̃(  , ) = 〈  ()〉 for the mean normalized correlation function of the back scattering phonon gas as a function of Knudsen number.

Fig. 7 .
Fig. 7. (a) Comparison of the normalized HCACF predicted by the multi-grey phonon gas models of SW Si (solid lines) with those from MD (translucent lines).The grey lines are for bulk Si while the coloured lines are for the pore geometries shown in the inset, where  = 54 nm,  = 26 nm and the pore radius is 5.4 nm.(b) Suppression function for the grey-phonon gas computed from the MC model with staggered pores, and the Matthiessen rule fit with  = 0.62.(c) The suppression function for the aligned pores, showing the Matthiessen fit with a crossover in scattering length from   = 0.92 in the diffuse regime to   = 0.17 in the ballistic regime.The inset plot in (b) shows the overlay of the suppression function for staggered and aligned pores.(d) Cumulative conductivity distribution over phonon mean-free-path normalized by   .(e) and (f) are scatter plots of phonon wave number vs suppression function for staggered and aligned pores respectively.The blue ovals indicate the long wavelength region where strong suppression is observed from the aligned pores, whereas the corresponding region in the staggered pore case is empty.