Introduction

In recent years, manipulating the phonon spectrum with two-dimensional (2D) phononic crystal (PnC) has attracted a great deal of research interest.1,2,3,4,5,6,7,8,9,10,11 Since different length scales are associated with phonons and electrons, thermal and electrical conductivity can be independently optimized in a periodic nanomesh of holey structure. Reproducible low-thermal conductivity, sufficient electrical properties and good mechanical strength make 2D Si PnCs as a promising candidate for thermoelectric applications, such as electrical power generation and on-chip thermal management for solid state devices.2 Although the impact of nano-scale periodic holes on the reduction of thermal conductivity is widely accepted, the underlying mechanism is still not clear. Generally, there are two contrasting pictures to describe phonon transport in PnCs, one is the particle-like incoherent phonon transport, and the other one is the wave-like coherent phonon transport. When the phonons lose their phase after successive scattering events, such as umklapp phonon–phonon scattering, phonon-impurity scattering and diffusive phonon-boundary scattering, these processes are often referred to as “incoherent” scattering mechanisms, in contrast to “coherent” scattering in which the phase is preserved and wave interference is occurred. Recently, the influence of coherent and incoherent scattering mechanisms on thermal conductivity of 2D PnCs is actively investigated in the field of nanoscale heat transport.12,13,14,15,16,17,18 Previous Boltzmann transport equation (BTE) calculations and Monte Carlo simulations considering the incoherent phonon transport largely overestimated the thermal conductivity of PnCs, compared to the experimentally measured results.3,4 Some groups ascribed the ultra-low-thermal conductivity of PnCs to coherent phonon interaction,1,5 which is analogous to that has been found in the superlattice structures.19,20

For coherent phonon transport, the periodicity of a hole array induces the folding of the Brillouin zone, modifies the phonon dispersion relation, and depresses phonon group velocities. It has been argued that coherent phonons affect thermal transport in periodic nanoporous films,1,3,4,5,12,13,14 and the range of occurrence of such phonon interference is a long-lasting question.13,14,15,16 In general, for coherent phonon modes associated with the secondary periodicity of the pores to affect thermal transport, the feature size of PnCs should be on a length scale comparable to wavelengths of the phonons that contribute to thermal transport.13 Zen et al. experimentally demonstrated the impact of the coherent scattering on thermal conduction of 2D PnCs at sub-Kelvin temperature, where phonon wavelengths were longer than the characteristic size of the structure.17 Maire et al. demonstrated that thermal conductivity was reduced in Si PnCs with an ordered array of holes as compared to that with randomly positioned holes below 10 K due to the presence of phonon interference, showed the temperature dependence of coherent effect and observed the transition from coherent to incoherent heat transport above 10 K.15 Furthermore, Wagner et al. fabricated PnCs with ordered and disordered lattices with equal filling fraction in free-standing Si membranes, and they found that the thermal conductivity of disordered PnCs was the same as that of ordered PnCs above room temperature.16 Using ultrafast pump-probe spectrum and Raman thermometry, they found that even in ordered PnCs, the phonon coherence only existed in the range of less than 0.4 THz at room temperature even for a hole wall roughness value as low as 1 nm. Thus significant modification of the thermal conductivity due to phonon coherence would only occur for very smooth surfaces/interfaces, where the frequency of the coherent phonon reaches the THz range, or at very low temperature, where the wavelength of the thermal phonon is remarkably increased. Considering the challenge in the fabrication of atomically smooth inner-wall and the strict periodicity, and the high-temperature environment for the application of thermoelectric energy generators, it is then conceivable that thermal conductivity of 2D PnCs in the incoherent regime plays critical role for practical thermoelectric devices. However, existing incoherent phonon models largely overestimate the thermal conductivity of 2D Si PnCs, hence the ultra-low thermal conductivity is still not well understood.

In this work, within the fully incoherent phonon picture, we develop a new model in which the phonon surface scattering due to bond order deficiency and the associated bond hardening in the surface skin of 2D PnCs is considered for the first time. We incorporate this strongly frequency-dependent phonon scattering process into the phonon BTE, and our theoretical calculations reproduce the ultra-low thermal conductivity of 2D Si PnCs reported in the recent experiments. Furthermore, the individual role of various phonon scatterings and their dominant frequency regimes in 2D Si PnCs are explored.

Results and discussions

Theoretical model and derivation

The termination of lattice periodicity in the surface normal direction of nanostructures will reduce the atomic coordination number (CN). Because of larger surface-to-volume ratio (SVR), a nanostructure has a high portion of under-coordinated atoms in the surface skin. Interatomic interaction and the changing fraction of the under-coordinated atoms are important factors for the different phonon transport behaviors of a nanomaterial from that of its bulk counterpart. Pauling and Goldschmidt indicated that, if the CN of an atom is reduced, the ionic and the metallic radius of the atom would shrink spontaneously.21,22 Therefore, the remaining bonds of the under-coordinated atoms become shorter and stronger.23 This effect deepens the potential well of trapping in the surface skin.24 Electrons in the relaxed region are more localized because of the depressed potential well of trapping, which lowers the conductivity in the surface region due to boundary scattering.25,26 Bond shortening and strengthening provide perturbation to the local potential, and this under-coordination effect should also have an important impact on the phonon transport dynamics. However, in previous model of phonon surface scattering, the change of interaction between under-coordinated atoms in the surface skin of nanostructures is overlooked. We present herein a novel phonon scattering process from the perspective of bond order imperfections in the surface skin on the basis of perturbation theory.

The spontaneous bond contraction and bond strength gain were formulated by Sun as follows:27,28

$$\left\{ {\begin{array}{*{20}{l}} {C_z = d_z{\mathrm{/}}d_0 = 2{\mathrm{/}}\left\{ {1 + {\rm exp}\left[ {\left( {12 - z} \right){\mathrm{/}}(8z)} \right]} \right\}} \hfill \\ {E_z = C_z^{ - m}E_{\rm b}} \hfill \end{array}} \right.$$
(1)

where C z is the bond-contraction coefficient, E z is the single-bond energy, z is the effective CN, d z is the bond length, Eb and d0 are the corresponding bulk values of single bond energy and bond length, respectively. m is a numerical exponent accounting for the nature of bond in a specific material, which is not freely adjustable. For Si the value of m is 4.88.29 From the perspective of dimensionality analysis, the force constant of bonds between under-coordinated atoms is given as

$$f_Z = {\rm d}^2u(r){\mathrm{/}}\left. {{\rm d}r^2} \right|_{r = {\rm d}_z} \propto E_z{\mathrm{/}}d_z^2$$
(2)

The bond order theory treats the outmost two atomic layers as the surface skin. The CN of atoms in the bulk phase, the outmost and the second outmost layer are 12, 4, and 6, respectively.27 According to Eqs. (1) and (2), the force constant between under-coordinated atoms in the surface skin is larger than that between atoms in the bulk. If there is an atom which has different force constant in crystal, the change of force constant δf results in the perturbation of Hamiltonian given as

$$H^{\prime} = \frac{1}{2}\delta f\mathop {\sum}\limits_{\bf{I}} {\left( {{\bf{u}}_{\bf{x}} - {\bf{u}}_{{\bf{x}} - {\bf{I}}}} \right)^2}$$
(3)

where I is a unit vector in the direction of a linkage, x is the lattice site of imperfection, and u x is the displacement of the xth atom from its equilibrium site. According to the second-order perturbation theory, the transition probability per unit time of an incident phonon of momentum \(\overrightarrow {\boldsymbol{k}}\) and energy ħω to a new state of momentum \(\overrightarrow {\boldsymbol{k}} ^\prime\) and energy ħω′ due to the perturbed Hamiltonian H′ is derived as30

$$p\left( {\overrightarrow {\boldsymbol{k}} ,\overrightarrow {\boldsymbol{k}} ^\prime } \right) = \frac{{2\left| {\overrightarrow {\boldsymbol{k}} \left| {H{\prime}} \right|\overrightarrow {\boldsymbol{k}} ^\prime } \right|^2}}{{\hbar ^2}}\left[ {\frac{{1 - {\rm cos}\left( {\omega {\prime} - \omega } \right)t}}{{\left( {\omega {\prime} - \omega } \right)^2t}}} \right]$$
(4)

Therefore the rate of phonon scattering is given as

$$\tau ^{ - 1} = {\int} {p\left( {\overrightarrow {\boldsymbol{k}} ,\overrightarrow {\boldsymbol{k}} ^\prime } \right)d\overrightarrow {\boldsymbol{k}} ^\prime }$$
(5)

From the perturbation theory of Klemens,31 the rate of phonon scattering by an isolated atom of different binding force in the bulk of crystal is derived as follows:

$$\tau (\omega )^{ - 1} = \frac{V}{{2\pi G^2v^3}}\left( {\frac{{\delta f}}{f}} \right)^2\omega ^4$$
(6)

If a number of imperfections are arranged at random in the bulk of crystal, the total scattering rate of phonons with frequency ω due to binding force imperfections, is proportional to the number of imperfections, which is given as

$$\tau (\omega )^{ - 1} = \frac{{VN_{\rm I}}}{{2\pi G^2v^3}}\left( {\frac{{\delta f}}{f}} \right)^2\omega ^4 = x\frac{V}{{2\pi Gv^3}}\left( {\frac{{\delta f}}{f}} \right)^2\omega ^4$$
(7)

where NI is the number of imperfections, f is the force constant, δf is the change of force constant, v is speed of sound, G is the number of atoms in the crystal, V is the volume of crystal, and x is the concentration of imperfections. The detailed derivation of Eqs. (6) and (7) is shown in Section 1 of Supporting Information.

The bond order imperfections are not randomly, but continuously arranged in the surface skin of nanostructures, therefore the scattering rate of phonons in Eq. (7) must be modified by two factors. Firstly, every bond links two under-coordinated atoms, so only half of under-coordinated atoms in the surface skin contribute to the perturbation of Hamiltonian, and the scattering rate in Eq. (7) must be multiplied by the factor of \(\frac{1}{2}\). Secondly, in the derivation of Eq. (6) in ref.31, every isolated atom of different binding force has six linkages, but the CN of atoms at the outmost atomic layer is 4 in our model, i.e., atoms in the outmost layer has only four linkages. Therefore, the scattering rate in Eq. (7) must be multiplied by the factor of \(\left( {\frac{4}{6}} \right)^2\). We then obtain the phonon scattering rate due to the bond order imperfections in the surface skin as follows:

$$\tau _{{\rm Bond}}^{ - 1} = \frac{1}{2}\left[ {x_1\left( {\frac{{f_4 - f_{\rm b}}}{{f_{\rm b}}}} \right)^2\left( {\frac{4}{6}} \right)^2 + x_2\left( {\frac{{f_6 - f_{\rm b}}}{{f_{\rm b}}}} \right)^2} \right]\frac{V}{{2\pi Gv^3}}\omega ^4$$
(8)

Here x1 = C4d0 · SVR (x2 = C6d0 · SVR) is the ratio between the number of atoms at the outmost (second outmost) layer and the total number of atoms. C4(C6) given by Eq. (1) is the bond-contraction coefficient of under-coordinated atoms at the outmost (second outmost) layer, d0 is the bond length in the bulk, and SVR is the surface-to-volume ratio. f4 (f6) given by Eq. (2) is the force constant between under-coordinated atoms at the outmost (second outmost) layer, and fb is the force constant between atoms in the bulk.

For silicon, d0 = 0.263 nm,32 v = 6400 m/s,33 V/G = a3/8, and the lattice constant a = 0.543 nm. From Eqs. (1), (2) and (8), the phonon scattering rate in silicon nanostructures due to bond order imperfections in the surface skin is

$$\tau _{{\rm Bond}}^{ - 1} = 1.86 \times 10^{ - 51} \cdot {\rm SVR} \cdot \omega ^4$$
(9)

where the unit of SVR is m−1. Because of the ω4 dependence, the phonon surface scattering by bond order imperfection is strongly frequency dependent. High-frequency phonons experience a much stronger surface scattering than low-frequency phonons. The traditional phonon-boundary scattering model is usually simply treated as frequency independent. By direct comparison with experiments, Wang et al. verified that this plausible treatment was not the best approximation for grain boundary scattering, and presented a new frequency-dependent model.34 Lee et al. experimentally found ballistic phonon transport prevailed in the cross-plane direction of holey silicon from 35 to 200 nm, and in their numerical model the lateral boundary scattering rate was treated by ω4, like Rayleigh scattering.9 Our theoretical work reveals the underlying mechanism of this frequency-dependent boundary scattering. The scattering rate in Eq. (9) is proportional to SVR, indicating that the importance of phonon-surface bond order imperfection scattering mechanism increases with decreasing the size of nanostructures. Actually, both Monte Carlo35 and molecular dynamics simulations36 showed that SVR is the universal gauge for thermal conductivity of nanowires.

We apply the phonon surface scattering to 2D Si PnCs with hexagonal lattice shown in Fig. 1. The temperature gradient is parallel to x-axis. According to the linearized phonon BTE within the relaxation time approximation, the thermal conductivity of phonons in branch λ in the x direction (direction of ∇T) is derived as37

$$\kappa _\lambda = \frac{V}{{(2\pi )^3}}{\int} {c_{{\rm ph}}v_x^2\tau d\overrightarrow {\boldsymbol{q}} }$$
(10)

where V is the volume of the sample, v x is the x component of the group-velocity vector, τ is the averaged phonon relaxation time between successive scattering events, \(\overrightarrow {\boldsymbol{q}}\) is the wave vector, and cph is the volumetric specific heat of each mode.

Fig. 1
figure 1

Schematic illustration of hexagonal-lattice Si PnCs. In our work, the thickness of 2D PnCs is fixed at 100 nm. The heat flux is along x direction. The neck size in our research is from 10 to 30 nm, and the pitch size is from 60 to 90 nm

Full size image

According to the Matthiessen’s rule (In Section 4 of Supporting Information, we evaluated the validity of Matthiessen’s rule), τ−1 = \(\tau _{\rm U}^{ - 1} + \tau _{\rm B}^{ - 1} + \tau _I^{ - 1} + \tau _{{\rm Bond}}^{ - 1}\), where \(\tau _{\rm U}^{ - 1}\), \(\tau _{\rm B}^{ - 1}\), \(\tau _{\rm I}^{ - 1}\) and \(\tau _{{\rm Bond}}^{ - 1}\) are the umklapp phonon–phonon scattering rate, traditional phonon-boundary scattering rate, phonon-impurity scattering rate and phonon-surface bond order imperfection scattering rate in Eq. (9), respectively. Slack and Galginaitis suggested the following form for the umklapp phonon–phonon scattering rate:38

$$\tau _{\rm U}^{ - 1}\left( \omega \right) = B_{\rm U}\omega ^2Te^{ - {\mathrm{\Theta }}/3T}$$
(11)

where Θ = ħωm/kB, is the Debye temperature, which is calculated by Brillouin zone-boundary frequency ωm.39

The traditional phonon-boundary scattering rate is given by

$$\tau _{\rm B}^{ - 1} = \frac{{1 - F}}{{1 + F}}\frac{{v(\omega )}}{{{\Lambda }_{\rm B}}}$$
(12)

where ΛB is the phonon-boundary mean-free path, F is the specularity parameter. Since in-plane scattering events by the PnC air holes and out-of-plane scattering by the PnC surfaces are not mutually exclusive processes but rather stochastic, \({\Lambda }_{\rm B}^{ - 1} = l^{ - 1} + {\Lambda }_{{\rm Pore}}^{ - 1}\),14 l is the thickness of PnC, ΛPore is the average distance between pore boundaries. In this work, fully diffusive scattering (specularity parameterF = 0) is assumed on the boundaries of 2D Si PnCs. In Section 3 of Supporting Information, we also investigate the influence of frequency-dependent diffusive/specular boundary scattering on the thermal conductivity of 2D Si PnCs. As shown in Fig. S4 in Supporting Information, when the surface roughness δ is larger than 0.7 nm, the results of frequency-dependent diffusive/partially specular boundary scattering are very close to that of fully diffusive boundary scattering. Therefore, it is reasonable to adopt the fully diffusive boundary scattering in this study as the reported surface roughness δ of 2D Si PnCs is usually larger than 1 nm.14,40

The phonon-impurity scattering rate is given by \(\tau _{\rm I}^{ - 1} = D\omega ^4\),41 the impurity scattering constant D for thin silicon films with varying doping concentrations is kept consistent with that of ref.12, i.e., 2.4 × 10−45, 4.0 × 10−45, and 7.0 × 10−44 s3 for low-doped, mid-doped, and high-doped holey silicon films.

The phonon-surface bond order imperfection scattering rate is given by Eq. (9), and SVR of hexagonal-lattice 2D PnC shown in Fig. 1 is derived as follows:

$${\rm SVR} = \left[ {\frac{{\pi \left( {p - n} \right)}}{{\frac{{\sqrt 3 }}{2}p^2 - \frac{\pi }{4}\left( {p - n} \right)^2}} + \frac{2}{l}} \right]$$
(13)

Here p, n, and l are the pitch size, neck size, and thickness of PnC, respectively.

Based on phonon BTE, the in-plane thermal conductivity (κ = κLA + 2κTA) of 2D Si PnCs can be derived as an integral, which is given by Eq. (10). The determination of all parameters, such as phonon dispersions and BU, Θ in Eq. (11) is in detail given in Section 2 of Supporting Information. The predicted thermal conductivities (κ) by Eq. (10) must be scaled by a volume correction function to account for the porosity ϕ in order to compare with the experimental data. The effective thermal conductivity is given by κEff = f(ϕ)κ.13,14 The volume correction function f(ϕ) = (1 − ϕ)/(1 + ϕ) was proposed by Hashin and Shtrikman,42 and agreed with the finite element method calculation by Jain et al.13 The porosity ϕ of hexagonal-lattice 2D PnCs shown in Fig. 1 is derived as follows:

$$\phi = \frac{{\frac{\pi }{4}\left( {p - n} \right)^2}}{{\frac{{\sqrt 3 }}{2}p^2}}$$
(14)

In-plane thermal conductivity of 2D Si PnCs

We investigate the in-plane thermal conductivity of 2D Si PnCs by incorporating our phonon-surface bond order imperfection scattering process to phonon BTE. Figure 2 compares the effective in-plane thermal conductivity of 2D Si PnCs with different neck size calculated by our model. The experimental data from ref.12 is also shown as benchmark. It is shown that the calculated results of our incoherent model agree well with the ultra-low thermal conductivity obtained in experiment. This verifies the concept proposed in ref.12 that the coherent phonon interference is disturbed because of the imperfections in the holey array and the resultant neck/pitch size error width (about 2 nm).

Fig. 2
figure 2

Effective in-plane thermal conductivity of 2D Si PnCs at different temperatures. The open circles present experimental values from ref.12, the solid lines present our calculated results. The red dotted line presents calculated results without surface bond order imperfection scattering mechanism. The doping level and neck size of red dot line are the same as those of red solid line. The pitch size is fixed at 60 nm

Full size image

In our phonon BTE calculations, excluding the phonon scattering due to surface bond order imperfections will result in large overestimation of thermal conductivity, which is demonstrated by the red dot line in Fig. 2. The doping level and neck size of red dot line are the same as those of red solid line. This demonstrates that the ignoring of phonon scattering by the surface bond order imperfections is the major source for the overestimation of thermal conductivity of 2D PnCs in other incoherent phonon transport models.

Individual role of various phonon scatterings in 2D Si PnCs

As shown in Fig. 2, in the case of fixed pitch size, κEff is dominated by neck size, and is nearly independent of the doping level. For example, κEff of low-doped 20 nm neck sized sample is lower than that of mid-doped 30 nm neck sized one, which indicates that the impurity scattering at these doping levels does not significantly affect the thermal conductivity of holey silicon. The umklapp peak in thermal conductivity of bulk silicon is absent, and smaller necks result in a weaker temperature dependence, which indicates that the umklapp phonon–phonon scattering is also not the dominant mechanism in governing thermal conductivity in nanoporous silicon films. The scattering rates of LA and TA phonons at room temperature as a function of phonon angular frequency are individually demonstrated in Fig. 3. It is clear that at room temperature the phonon surface scattering by bond order imperfections dominates over other scatterings when ω is higher than 35.7 Trad/s for LA phonons and 25.9 Trad/s for TA phonons. Below 35.7 Trad/s for LA phonons (or 25.9 Trad/s for TA phonons), the traditional phonon-boundary scattering is the dominant mechanism. Not only the underlying mechanism but also the frequency dependence of phonon-surface bond order imperfection scattering and traditional phonon-boundary scattering are quite different, therefore both of them must be considered in the phonon transport of nanostructures. The suppression of high-frequency phonons by surface bond order imperfections scattering in company with the suppression of low-frequency phonons by normal boundary scattering result in the remarkable reduction of thermal conductivity of 2D Si PnCs. According to the phonon dispersions of LA and TA branches used in our calculation and Bose–Einstein distribution of phonons, there are 79.8% LA phonons with ω > 35.7 Trad/s, and 64.8% TA phonons with ω > 25.9 Trad/s at room temperature.

Fig. 3
figure 3

Various scattering rates of LA phonons (a) and TA phonons (b) as functions of phonon angular frequency at room temperature. The pitch size of holey silicon is 60 nm, and the neck size is 20 nm

Full size image

The severe scatterings of high-frequency phonons by surface bond order imperfections evidently change the normalized accumulative distribution of thermal conductivity, as shown in Fig. 4. Here, Φ(ω) is the relative contribution to the thermal conductivity for phonons with angular frequency from 0 to ω. The solid line shows the results based on our new model, and the dotted line presents the results from the model without the phonon scattering by surface bond order imperfections. Because high-frequency phonons suffer much stronger scattering by surface bond order imperfections than low-frequency phonons, their contribution to thermal conductivity is remarkably suppressed. For example, with consideration of surface bond order imperfections scattering, the phonons of LA branch with frequency ω > 60 Trad/s have about 24% contribution to thermal conductivity at room temperature, but phonons in the same frequency regime have 46% contribution to thermal conductivity if this surface scattering process is not considered. The overestimation of thermal conductivity in other incoherent models stems from the lack of severe scattering of high-frequency phonons by surface bond order imperfections.

Fig. 4
figure 4

Normalized accumulative distribution of thermal conductivity of LA phonons (a) and TA phonons (b) as a function of phonon angular frequency for holey silicon at room temperature. The pitch size of holey silicon is 60 nm, and the neck size is 20 nm

Full size image

Very recently, Lee et al. found experimentally the identical thermal conductivities for periodic and aperiodic nanomeshs (NMs) of the same porosity and average pitch size.43 Other works such as refs.15,16 demonstrated the similar observation. It is shown in Fig. 3 that phonon-surface bond order imperfection scattering is the dominant process for thermal phonons. According to our theory, the scattering rate of this process is proportional to SVR, given in Eq. (9). In the case of same feature size and porosity, SVR of ordered and disordered NMs is equal, therefore the thermal conductivities are identical.

Neck-size and pitch-size dependence of thermal conductivity

Furthermore, ignoring the phonon scattering by surface bond order imperfections also results in discrepancy from experimentally observed neck-size dependence in thermal conductivity. Figure 5 shows the dependence of thermal conductivity of holey silicon on neck size and pitch size at room temperature. Both increase of pitch size and decrease of neck size lead to decrease of κEff, but the influence of neck size on κEff is stronger than that of pitch size. A correlation could be captured using a power law fit, κEff ∝ nm and κEff ∝ pm. The fitting constant m is found to be 1.9 for neck size and −1.2 for pitch size. The solid lines in Fig. 5 are the best-fitting results. By experimental measurements, Lim et al. found power law relationship between κEff and neck size (fixed pitch size), and the neck size-dependent scaling factor (m) was found to be ~2,12 which is very close to our calculation result. In our calculation, if we ignore the phonon scattering by surface bond order imperfections, the neck size-dependent scaling factor m is close to 1. Therefore, our theoretical work indicates that the phonon-surface bond order imperfection scattering plays the dominant role in obtaining quantitative agreement with the experimental observations of the neck-size and pitch-size dependence of the thermal conductivity.

Fig. 5
figure 5

Dependence of effective thermal conductivity of holey silicon on neck size (a) and pitch size (b) at room temperature. The data dots are our calculated results, the red line is the best-fitting curve κEff ∝ n1.9, and the blue curve is the best-fitting relationship κEff ∝ p−1.2

Full size image

Effects of roughness of the hole wall on thermal conductivity

It is worth mentioning that achieving ultra-low-thermal conductivity in the incoherent phonon transport limit is critically important for further reducing thermal conductivity of 2D Si PnCs. For representative nanostructures, such as SiNWs, introduction of surface roughness is a commonly adopted method to reduce thermal conductivity.35,44,45,46,47,48,49,50,51 However, from the viewpoint of coherent models, this method (introduction of surface roughness) is believed to be a “negative” effect to reduce thermal conductivity of 2D PnCs, because it will destroy the coherence of phonons. In contrast, as demonstrated in our present work, phonon coherent transport is not the necessary condition to achieve ultra-low thermal conductivity in 2D PnCs. Therefore, surface engineering becomes a promising choice to further reduce the thermal conductivity of 2D PnCs. In our model, it is clear that SVR is the key factor which affects phonon transport and thermal conductivity in holey silicon. As demonstrated in Fig. 3, surface bond order imperfection scattering is a filtering mechanism for high-frequency phonons, and it is stronger with larger SVR according to Eq. (9). Hence, our theoretical model indicates that roughening the hole wall to increase SVR can further reduce the thermal conductivity of 2D Si PnCs. Next we will explore the effects of roughness of the hole wall on the thermal conductivity of 2D Si PnCs.

It is general to assumed that the autocorrelation function (ACF) of a random rough surface follows an exponential curve or Gaussian curve (exponential or Gaussian surfaces, for brevity),52 which is given by

$$C\left( {\vec r} \right) = \left\{ {\begin{array}{*{20}{l}} {\delta ^2{\rm e}^{ - \vec r/\xi }} \hfill & {{\mathrm{exponential}}} \hfill \\ {\delta ^2{\rm e}^{ - \vec r^2/\xi ^2}} \hfill & {{\mathrm{Gaussian}}} \hfill \end{array}} \right.$$
(15)

where δ denotes the root-mean-square (rms) value of roughness, ξ is the correlation length which is related to the lateral length scale of the roughness. By convolution theorem and inverse Fourier transform,53 we can generate real-space rough hole surface profile, and then the surface area of holey silicon and SVR can be calculated. ξ is selected to be 6 nm in our calculations. As shown in Fig. 6a, increasing the roughness of hole wall effectively decreases the thermal conductivity of holey silicon. In the case of same roughness, the reduction of thermal conductivity by exponential rough surface is more significant than that by Gaussian surface, because the exponential surface has considerably more small-scale roughness shown in Fig. 6b, which results in larger SVR than that of Gaussian surface.

Fig. 6
figure 6

Influence of roughness on κEff and SVR of holey silicon nanostructures at room temperature (a). Rough surface profiles of hole wall generated from Exponential and Gaussian ACF (b). R is the radius of hole. The pitch size of holey silicon is 60 nm, and the neck size is 30 nm

Full size image

In conclusion, we have presented a novel phonon scattering process, which stems from bond order imperfections in the surface skin of nanostructures. We have incorporated this strongly frequency-dependent scattering process to phonon BTE, and reproduced the ultra-low experimental thermal conductivity of holey silicon nanostructures. According to our analysis, the suppression of high-frequency phonons by surface bond order imperfections scattering in company with the suppression of low-frequency phonons by normal boundary scattering result in the remarkable reduction of thermal conductivity. We have found that the filtering mechanism for high-frequency phonons is proportional to SVR. It is demonstrated that increasing SVR by roughening the hole wall is an effective approach to further reduce the thermal conductivity of 2D Si PnCs. Our model is helpful for understanding the limitation of surface on phonon transport, as well as for devising novel methods to manipulate thermal transport in nanostructures by surface engineering. This model is general for PnCs and can be extended to nanostructures of various materials.

Data availability

The data that support the findings of this study and the code for the strategy proposed in this study are available from the corresponding author (G. F. Xie) upon reasonable request.