Influence of point defects and multiscale pores on the different phonon transport regimes

A common strategy to tailor the thermal conductivity of a material is to introduce structural features that modulate phonon scattering, such as atomic-scale defects and nano- and macro-sized pores. However, particle-like and wave-like phonon transport and scattering during a crossover in thermal transport regimes is not well understood. Here, we perform a rigorous quantitative comparison of the thermal conductivity obtained from molecular dynamics simulations and phonon Boltzmann transport equations, taking graphene as an example. We observe a generally increasing trend in thermal conductivity when the pore size increases from point defect to nanopore, due to a transition from Rayleigh scattering to geometric scattering and reduced boundary density. The thermal conductivity further converges to the diffusive limit for macropores because of the dominant effect of phonon-phonon scattering over phonon-boundary scattering. Moreover, we identify a critical interpore distance for the crossover from dependent to independent phonon-pore scattering and a critical pore size for the crossover from point defect scattering to boundary scattering. This work provides a comprehensive understanding of phonon transport in materials containing defects and pores. Structural features control the thermal conductivity of a material by modulating phonon scattering. Here, simulations and theory reveal the effect that atomic-scale defects and pores have on the crossover of thermal transport regimes in graphene.

lectromagnetic wave scattering by the particle system has been widely studied owing to its great importance to science and engineering [1][2][3][4] . For an individual particle, the scattering includes Rayleigh, Mie, and geometric scattering regimes 1,5 . For multiple particle ensembles, the interactions of the particles can further manifest, such as multiple or dependent scattering 6 . In comparison, lattice waves, i.e., thermal phonons are broadband in the spectrum and have strong anharmonic interaction with others 7 . The transport and scattering processes of phonon with imperfections are probably equally complicated, if not more. By virtue of recent advances in material synthesis 8 , it is now possible to introduce pores of different scales ranging from atomic-scale vacancies to macropores in a controllable manner [9][10][11] . On the practical level, further understanding of phonon scattering with pores of different scales is particularly important for modulating the thermal conductivity of materials for different applications, like thermoelectric energy conversion [12][13][14][15] .
Phonons can behave like particles or waves depending on circumstances 16 . In the past decades, many efforts have been made to establish the theory of phonon scattering with imperfections based on the particle-like (incoherent) or wave-like (coherent) characteristics of phonons [17][18][19][20][21][22][23][24][25] . By analogy with the electromagnetic wave 6,26 , particle-based phonon radiative transport equation has been proposed 21 . Mie scattering theory has been developed 18,23 to consider the wave-like dependent scattering effects 18,22 in phonon transport. These theoretical models are derived based on a single (usually long) wavelength and are difficult to consider the collective behaviors of broad-spectrum phonons in real materials, e.g., thermal conductivity. On the other hand, simulation methods have been developed to explore the thermal transport mechanisms in porous materials by examining the impact of porosity 27 , pore shape 28 , and pore arrangement 29 on thermal conductivity. Nevertheless, previous works usually focus on a single pore scale: atomic defect, nanopore, and macropore, which is usually dealt with atomistic simulation methods 30,31 , phonon Boltzmann transport equation (BTE) 32 , and heat diffusion equation 33 , respectively. The crossovers of thermal transport for different pore scales have not been clearly identified. Especially, the interplay of particle-like and wave-like phonon transport and scattering in the crossovers is not uncovered.
In this paper, we quantitatively study the crossover of different phonon transport regimes with pores of different scales ranging from point defect to macropore employing graphene as the model system owing to its scientific significance [34][35][36][37] . Due to the broadband nature of phonons, we focus more on the collective effects of phonon transport on thermal conductivity, instead of single phonon mode analysis 22 . The key idea of distinguishing different phonon transport regimes is to decompose the wave-like and particle-like phonon transport by comparing the thermal conductivity computed from molecular dynamics (MD) simulation and phonon BTE. Unlike electromagnetic waves, which can be described by the wave equation, i.e., Maxwell equation 1 and the transport equation, i.e., radiative transport equation 38 , a major challenge to studying phonon transport is the lack of a general governing equation for wave-like phonon transport. So far, the best tool that can include the wave nature of phonon is the MDs simulation 39 . The phonon transport excluding the wave effect can be well described by phonon BTE 40 . A quantitative agreement of MD and phonon BTE is recently enabled by an indepth understanding of MDs simulation with phonon picture 41 . It has been shown that carefully-implemented MD and wellconverged mode-resolved phonon BTE yield quantitatively identical thermal conductivity values for silicon thin films 41 . This provides a solid ground to isolate the wave-like and particle-like phonon transport 42 and further identify the crossover of different phonon transport regimes in porous materials.
In order to study the phonon-pore scattering at different scales and eliminate the phonon ballistic transport effect, the system must be "bulk", i.e., the system size must be much larger than the pores and the phonon mean free path. To achieve that, we study a periodic simulation domain, which is an infinitely large system. Accordingly, we employ homogenous nonequilibrium molecular dynamics (HNEMD) simulation 43 , which has been recently implemented in an efficient GPU-based MD code 44 and has been shown more efficient than the equilibrium MDs simulation method 45 . This allows us to study a domain up to 103 nm by 102 nm with reasonable computational cost. The simulation setup of HNEMD for the porous graphene is shown in Fig. 1a. We consider the thermal conductivity along the y-direction. Periodic boundary conditions 46 are applied along the x and y-directions. L and W denote the domain length and width, respectively. D denotes the diameter of the pore. The number of absent atoms of the pore is denoted as N. The temperature is set as 300 K. Correspondingly, the periodic simulation domain is considered in phonon BTE calculation. Figure 1b shows the simulation setup of phonon BTE calculation for porous graphene of the same geometry as in HNEMD. Periodic boundary conditions are applied at the top, bottom, left, and right boundaries of the simulation domain. The boundary of the pore is set as a diffusely reflecting boundary 47 . In order to make a fair quantitative comparison, the intrinsic phonon properties of graphene involved in HNEMD and BTE are ensured the same (Method section). Through a rigorous comparison between the thermal conductivity computed by HNEMD and phonon BTE, we identify the crossover from dependent to independent scattering and the crossover from point defect scattering to boundary scattering in porous graphene. Furthermore, the impact of pore size crossing different scales on the thermal conductivity of porous graphene in the independent and dependent scattering regimes is investigated. Finally, the thermal transport regimes of point defect, nanopore, and macropore are distinguished.

Results and discussion
Crossover of dependent and independent phonon-pore scattering. We first study the crossover from dependent to  Fig. 1 The schematic of isolating wave-like and particle-like thermal transport by comparing HNEMD and phonon BTE simulation. a The setup of HNEMD simulation. b The setup of phonon BTE calculation. Periodic boundary conditions are applied to x-directional and y-directional boundaries in both HNEMD and BTE. In BTE, the boundary of the pore is set as a diffusely reflecting boundary. D denotes the diameter of the pore. The number of absent atoms of the pore is denoted as N. L and W denote the domain length and width, respectively. The temperature is set as 300 K. The thermal conductivity is along the y-direction.
independent phonon-pore scattering in porous graphene for different scales of pores. For electromagnetic waves, according to the volume fraction and the particle size relative to the wavelength, the scattering in discrete random media can be divided into independent and dependent scattering regimes 26,48 . In thermal transport, dependent and independent scattering regimes also exist for phonon-pore scattering but the crossover is never examined. As shown in Fig. 2a, when the interpore distance δ is comparable to the wavelength of phonon λ, i.e., δ~λ, dependent phonon-pore scattering will prevail, which is originated from interference between waves scattered from nearby pores 22,49 . If the interpore distance δ is much larger than the phonon wavelength λ, i.e., δ>>λ, the interference effects could be neglected and phonon-pore scattering is independent 22 , as presented in Fig. 2b.
To observe the crossover from dependent to independent phonon scattering, we consider one pore located at the center of the simulation domain with increasing interpore distance at a fixed pore size, as shown in Fig. 2c. Periodic boundary conditions are applied to the boundaries of the domain (denoted by the dashed line). The interpore distance increases from 3 nm to 80 nm. Three different diameters of the pore are considered: 0.8 nm, 5 nm, and 20 nm, which span the range of dominant phonon wavelength of graphene at room temperature, as shown in Supplementary Fig. 1 in the Supplementary Information. As will be shown later, the pore size crosses from the scale of point defect to nanopore. Figure 2d presents the thermal conductivities of porous graphene at different interpore distances computed from HNEMD simulation and phonon BTE calculation for three different diameters. It is seen that thermal conductivities computed from HNEMD (κ MD ) and phonon BTE (κ BTE ) both increase with increasing interpore distance because of the decreasing porosity. We then quantify the relative difference between κ MD and κ BTE as η MD-BTE = (κ MD -κ BTE )/κ MD × 100%, which is shown in Fig. 2e. It is found that for all three different diameters, the absolute value of the relative difference |η MD-BTE | generally decreases with increasing interpore distance and converges at the interpore distance of around 10 nm.
The deviations between κ MD and κ BTE come from the wave effects that are not considered in BTE, which consist of two parts. One is the wave-like phonon scattering by an individual pore, which manifests at a small pore size. Another part is the wave interference of scattered phonons between nearby pores, which increases with decreasing interpore distance. Since the diameter of the pore is fixed, the deviation between κ MD and κ BTE from the wave-like phonon scattering by an individual pore can be regarded as unchanged as the interpore distance increases. Thereafter, the change of deviation η MD-BTE with interpore distance mainly originated from the modification of the wave interference effect. When the interpore distance is small, strong wave interference occurs and dependent scattering dominates. With increasing interpore distance, the wave interference is weakened and phonon-pore scattering transits from dependent to independent. Therefore, the change of deviation η MD-BTE with increasing interpore distance can indicate a crossover from dependent scattering to independent scattering. As indicated in Fig. 2e, when the interpore distance is roughly larger than 10 nm, the phonon scattering with pores can be regarded as independent for graphene at room temperature. According to Supplementary Note 1 in the Supplementary Information, the dominant phonon wavelength of graphene at room temperature is~1 nm. Therefore, we may conclude that when the interpore distance is larger than ten times the dominant phonon wavelength, the dependent scattering between pores will have a negligible effect on the thermal conductivity. Although there could be some phonons of long wavelengths that still experience-dependent scattering, the contribution of these phonons to the total thermal transport is small, and the overall phonon-pore scattering can be still regarded as independent.
Crossover of point defect scattering and boundary scattering. We then consider the phonon scattering with an individual pore. For electromagnetic waves, based on the ratio between wavelength and particle size, the scattering cross section can be divided into three regimes where the wavelength is larger, roughly equal, and smaller than the particle size 49,50 . In the first regime, the cross-section obeys Rayleigh scattering law as the fourth power to frequency (~ω 4 ) 50 . In the third regime, the cross-section follows the frequency-independent geometric scattering 50 . In the second regime, the scattering cross section is in the Mie scattering regime, which can have oscillations with wavelength 50 . Similarly, Interpore distance (nm) Interpore distance (nm) The crossover from dependent to independent phonon-pore scattering. a The schematic of dependent phonon-pore scattering. Dependent phonon-pore scattering refers to the situation when the interpore distance is comparable to the phonon wavelength (δ~λ) that the interference among scattered waves from nearby pores will prevail. b The schematic of independent phonon-pore scattering. Independent phonon-pore scattering refers to the situation when the pores are sufficiently distant (δ >> λ) that the interference effect can be safely neglected. c The simulation setup for studying the crossover from dependent to independent phonon-pore scattering. There is one pore located at the center of the simulation domain with increasing interpore distance. The pore size is fixed. Periodic boundary conditions are applied to the boundaries of the domain. d The thermal conductivities of porous graphene at different interpore distances computed by HNEMD and phonon BTE. e The relative difference between HNEMD and BTE computed thermal conductivities with respect to the interpore distance. The light blue and light yellow regions denote the dependent and independent regimes, respectively.
there exists the crossover of Rayleigh scattering to geometric scattering for phonon scattering with imperfections, but it is less explored.
Owing to the broadband nature of phonons, it is necessary to consider the phonons with different wavelengths. When the entire phonon wavelength λ is larger than the pore size, the pore can be regarded as a point defect, and the phonon scattering with the pore generally obeys the law of Rayleigh scattering [51][52][53] , as shown in Fig. 3a. In contrast, when the entire phonon wavelength λ is much shorter than the pore size, phonon scattering with pores generally follows the geometric scattering, which can be described as frequency-independent boundary scattering with a nanopore, as shown in Fig. 3c. When the pore size is approximately comparable to some phonon wavelengths, phonons with longer wavelengths than the pore size would experience Rayleigh-type scattering, whereas phonons with shorter wavelengths are scattered geometrically, as shown in Fig. 3b. In all the regimes, there are always phonons experiencing Mie scattering. It is difficult to distinguish the Rayleigh or Mie scattering only based on the thermal conductivity difference because the wave effects are both significant in the two regimes. Therefore, we distinguish the overall phonon scattering with pores as point defect scattering or boundary scattering of nanopore rather than Rayleigh, Mie, or geometric scattering regimes.
To observe the crossover from point defect scattering to boundary scattering, we consider one pore located at the center of the simulation domain with increasing pore size at a fixed interpore distance, as shown in Fig. 3d. Periodic boundary conditions are applied to the boundaries of the domain (denoted by the dash line). The pore size varies from 0.17 nm (N = 1) to 10.22 nm. To exclude the interference effect, the interpore distance is set as 20 nm and 30 nm based on the results presented in Fig. 2e. κ MD and κ BTE of porous graphene with increasing pore size for the interpore distance of 20 nm and 30 nm are shown in Fig. 3e, which both decrease with increasing pore size because of the increasing porosity. The relative difference between κ MD and κ BTE is shown in Fig. 3f. It can be seen that for both δ of 20 nm and 30 nm, with increasing pore size, the relative difference η MD-BTE generally increases from negative to zero and converges at 4% at the diameter of 5 nm.
Since the interference effect can be neglected at an interpore distance larger than 10 nm, the deviation between κ MD and κ BTE only results from the phonon scattering by an individual pore. MD includes both the wave-like scattering of long-wavelength phonons and particle-like scattering of short-wavelength phonons. On the contrary, phonon BTE takes phonons as particles, in which the wave-like phonon scattering is ignored and only particle-like geometric scattering is considered 29 . Since the Rayleigh scattering cross section is much larger than that of geometric scattering 20 , BTE overestimates the thermal conductivity for small pores owing to the neglect of the wave-like phonon-defect scattering. With increasing pore size, the particlelike phonon-boundary scattering dominates over the wave-like phonon-defect scattering, and the difference between κ MD and κ BTE decreases. Therefore, the variation of η MD-BTE with increasing pore size can indicate a crossover from the regime dominated by point defect scattering to boundary scattering. As shown in Fig. 3f, phonon-pore scattering is dominated by boundary scattering when the pore size is larger than 5 nm, which is around five times the dominant phonon wavelength, as shown in Supplementary Fig. 1 in the Supplementary Information.
It is noted that the difference between κ MD and κ BTE in the boundary scattering regime is slightly larger than zero. We believe that the small deviation might stem from the assumption that the phonon scattering with a pore in MD could be equivalent to phonon scattering by a diffusive boundary in BTE. Theoretically, there is probably no rigorous correspondence of phonon scattering with pores described by the MD simulation. Since the concept of "phonon-boundary scattering" is already based on  Fig. 3 The crossover from point defect scattering to boundary scattering. a The schematic of point defect scattering. When the entire phonon wavelength λ is larger than the pore size, the pore can be regarded as a point defect, and the phonon scattering follows Rayleigh scattering. b The schematic of phonon scattering between point defect scattering and boundary scattering. When the pore size is comparable to some phonon wavelengths, phonons with longer wavelengths than the pore size would experience Rayleigh scattering, whereas phonons with shorter wavelengths are scattered geometrically. c The schematic of boundary scattering. When the entire phonon wavelength λ is much shorter than the pore size, phonon scattering with pores follows geometric scattering, which can be described as phonon-boundary scattering. d The simulation setup for studying the crossover from point defect scattering to boundary scattering. There is one pore located at the center of the simulation domain with increasing pore size. The interpore distance is fixed. Periodic boundary conditions are applied to the boundaries of the domain. e Thermal conductivities of porous graphene with different pore sizes at the interpore distance of 20 nm and 30 nm were obtained by HNEMD and phonon BTE. f The relative difference between HNEMD and BTE computed thermal conductivities for different pore sizes. The light blue and light orange regions denote the thermal transport regime dominated by phonon-defect scattering and boundary scattering, respectively. The error bar denotes the standard deviation.
the particle transport picture and phonons are broad-spectrum, it is never possible to set boundary conditions at the pore edge in BTE that is completely equivalent to MD. Nevertheless, the general tendency of deviation η MD-BTE could still provide evidence for the crossover from point defect scattering to boundary scattering.
Impact of pore size on the thermal conductivity. We further investigate the impact of pore size crossing different scales on the thermal conductivity of porous graphene. The simulation domain size and the porosity are fixed. Porous graphene with periodic and aperiodic pore arrangements is considered, of which the schematic atomic configurations are shown in Fig. 4a, b, respectively. The aperiodic configuration is generated by randomly shifting the centers of pores from their periodic positions 42 . Moreover, the independent and dependent scattering regimes are considered, of which the simulation setup is shown in Fig. 4c, d, respectively. For the independent scattering regime, the simulation domain is set as 72.3 nm × 71.6 nm, and the porosity is fixed as 0.02%. The pore diameter increases from 0.17 nm (N = 1) to 1.10 nm (N = 36). With increasing pore size, the corresponding interpore distance increases, which is always larger than 10 nm (ensuring independent scattering). For the dependent scattering regime, the simulation domain is set as 26.6 nm × 23 nm, and the porosity is fixed at 2.8%. The diameter increases from 0.17 (N = 1) nm to 4.64 nm (N = 648), which crosses from point defect to nanopore regime as indicated in Fig. 3f. The interpore distance is always smaller than 10 nm, ensuring the dependent scattering. For both two scenarios, periodic boundary conditions are applied to the boundaries of the domain (denoted by the dashed line). Figure 4e, f presents thermal conductivities of periodic and aperiodic porous graphene with increasing pore sizes computed by HNEMD and phonon BTE for the independent and dependent scattering regimes, respectively. The thermal conductivities of periodic and aperiodic porous graphene are denoted as "p" and "ap" in Fig. 4e, f. The BTE computed thermal conductivities of periodic and aperiodic porous graphene are almost exactly the same, and thus only the results of periodic porous graphene are presented, which are denoted as "BTE". As shown in Fig. 4e, the thermal conductivities of periodic and aperiodic porous graphene are quite close, which indicates that phonon scattering by each pore is independent and not affected by other pores. BTE generally overestimates the thermal conductivity due to the neglect of wave-like phonon-defect scattering as mentioned before. κ BTE monotonically increases with pore size because of the reduced boundary density 54 . κ MD generally first decreases then increases with the pore size as denoted by the magenta dash line, which agrees with the tendency of phonon scattering with nanoparticles 55,56 . Mingo et al. predicted a minimum of thermal conductivity with an increase in nanoparticle size for the alloy system with nanoparticles 57 . In their later work 58 , they discovered that for graphene with isotope clusters at cryogenic temperatures, the thermal conductivity decreases with the cluster size. They claim that for larger clusters, the thermal conductivity starts to increase 58 . At a small pore size, the scattering cross section increases with the pore size following the Rayleigh scattering mechanism, and thus the thermal conductivity decrease with increasing pore size 59 . On the other hand, at a large pore size, the scattering cross-section is proportional to the geometric cross-section, and the thermal conductivity increase with pore size 59 . The total scattering cross section is the interpolation of Rayleigh and geometric scattering cross sections 20 . Accordingly, it is expected that thermal conductivity first decreases and then increases with the pore size.
It is noted that the κ MD of D = 0.17 nm (N = 1) and 0.45 nm (N = 6) are far from the trend. This is attributed to the fact that the atomic configuration of the pore is not circular and thus the changes in scattering cross section do not rigorously correspond to the theoretical predictions 59 . To clarify, the atomic configurations of the pore of N = 1, 2, 3, 4, 6, 9 are presented in Fig. 4e. The thermal conductivity is calculated along the y-direction. It is seen that the cross section perpendicular to the y-direction of a pore of N = 1 is the same as that of N = 2, which is 2 ffiffi ffi 3 p times the lattice constant. Therefore, the scattering cross-section of a pore of N = 1 is close to that of N = 2. Whereas, the number of pores for N = 1 is larger than that for N = 2. Therefore, the total cross sections of N = 1 are much larger than that of N = 2 and the thermal conductivity of N = 1 is lower than that of N = 2. This phenomenon is also explained by the difference in undercoordinated atoms between single vacancy and double vacancy 60 . Similarly, the cross-section perpendicular to the y-direction of N = 6 is close to that of N = 4, while the number of pores for N = 6 is smaller than N = 4. Thus, the total cross sections of N = 6 are much smaller than that of N = 4 and the thermal conductivity of N = 6 is higher than that of N = 4 and even larger than that of N = 9.
Then we consider the dependent scattering regime as shown in Fig. 4f. The thermal conductivities of aperiodic porous graphene are much lower than that of periodic counterparts, indicating the existence of phonon localization induced by dependent scattering 61,62 , which is confirmed by the spectral analysis of phonon density of states 39 and participation ratio 39 in Supplementary Note 4 in Supplementary Information. BTE overestimates the thermal conductivity of aperiodic structures and underestimates the thermal conductivity of periodic structures for smaller pores. Then, the deviations decrease with increasing pore size due to the crossover from point defect scattering to boundary scattering, as indicated in Fig. 3f. κ MD of periodic structures first decreases and then increases with the pore size, similar to the trend for the independent scenario shown in Fig. 4e. Whereas, κ MD of aperiodic structures monotonically increases with pore size, which indicates the dominant role of particle-like transport after the cancellation of wave-like Rayleigh scattering and phonon localization 42 .
Thermal transport regimes of point defect, nanopore, and macropore. Lastly, we distinguish the thermal transport regimes of point defect, nanopore, and macropore. For nanoporous materials, the thermal transport is ballistic, which can be described by the phonon BTE when the wave effects can be ignored. The thermal conductivity of nanoporous materials is dependent on the pore size as indicated in Fig. 4. As for macroporous materials, the thermal transport is diffusive, which can be described by Fourier's law. The thermal conductivity of macroporous materials is independent of the pore size. The comparison between the thermal conductivity predicted by phonon BTE and Fourier's law can be used to identify the crossover from nanopore to macropore.
Taking porosity of 2.8% as an example, we compute the thermal conductivities of porous graphene with different pore sizes ranging from several nanometers to tens of micrometers by phonon BTE and Fourier's law. The results are shown in Fig. 5. As the pore size increases, the BTE computed thermal conductivity first increases due to the reduced boundary density and gradually converges owing to the dominance of phononphonon scattering over the phonon-boundary scattering. The converged value for BTE corresponds to the thermal conductivity value predicted by Fourier's law, which occurs at the pore sizes 1 μm. The corresponding interpore distance is~5 μm. For graphene at room temperature,~85% of the heat is carried by the phonons of the mean free paths shorter than 1 μm 42 . Therefore, it is reasonable that porous graphene with pore size larger than 1 μm can be regarded as macroporous materials. Based on the two critical pore sizes of 5 nm (determined before) and 1 μm, the thermal transport regimes of point defect, nanopore, and macropore of graphene can be distinguished, which are indicated Interpore distance (nm) Point defect Macropore Nanopore in Fig. 5. For the point defect regime, the thermal conductivities computed by MD simulations have been presented in Fig. 4f, and the same data are replotted in Fig. 5.
As seen in Fig. 2, the interpore distance has an important effect on the thermal conductivity of porous graphene. The interpore distance has been used to account for the phonon size effect induced by boundary scattering in the kinetic relationship 14,63,64 κ ¼ 1 2 FðφÞ R CvΛ eff dΛ. C is the phononspecific heat, v is the phonon group velocity and F(φ) is the correction factor for the porosity φ. Λ eff is the effective mean free path, accounting for the phonon size effect induced by boundary scattering, which can be modified from the bulk phonon mean free path (MFP) Λ bulk based on Matthiessen's rule as Λ eff = (1/Λ bulk + 1/Λ pore ) −1 . Λ pore is the characteristic length of the porous structure 14,63,65,66 . We use the kinetic model to predict the thermal conductivities of porous graphene by taking Λ pore as interpore distance (more details can be found in Supplementary Note 5 in Supplementary Information), which is shown as the solid line in Fig. 5. The predicted thermal conductivity by the kinetic model agrees well with the BTE calculation, indicating that the interpore distance can truly capture the size effect induced by phonon-boundary scattering.

Conclusion
In summary, we investigate the phonon transport in porous graphene with pores of different scales ranging from point defects to macropores. Through a rigorous quantitative comparison between the thermal conductivity computed by HNEMD simulation and phonon BTE, we identify the critical interpore distance for the crossover from dependent to independent phonon scattering and the critical pore size for the crossover from point defect scattering to boundary scattering, which is around ten and five times the dominant phonon wavelength, respectively. We find that the thermal conductivity of porous graphene generally increases and converges to the diffusive limit as the pore increases from point defect to macropores. The fundamental mechanisms include the transition from Rayleigh to geometric scattering in the point defect regime, the reduced boundary density in the nanopore regime, and the dominant effect of phonon-phonon scattering over the phononboundary scattering in the macropore regime. Our work facilitates a comprehensive understanding of phonon scattering with pores and thermal transport in defected materials, which can provide theoretical guidance for the engineering of multiscale defects to manipulate the thermal conductivity of materials.

Methods
MDs simulation. In HNEMD, the simulated system is in a non-equilibrium state, which is generated by the external driving force homogeneously applied to the system rather than two thermal baths at different temperatures 67 . As such, the boundary scattering in the transport direction is absent in the HNEMD simulation 46 . The details of the derivation of the thermal conductivity expression in the HNEMD method can be found in ref. 45 . Here, we only present the final expression of average running thermal conductivity κ ave (t), which is given as 45 J q is the heat current, T is the temperature, V is the system volume, F e is the driving force parameter and 〈〉 ne denotes the nonequilibrium ensemble average. The integral represents a running average of thermal conductivity over a simulation up to time t. The simulation setup of HNEMD is shown in Fig. 1a. We consider the thermal conductivity along the y-direction. Periodic boundary conditions 46 are applied along the x and y-directions. All the HNEMD simulations in this work are performed with GPUMD package 44 . To model the thermal transport, Tersoff potential is applied with the parameters from ref. 68 . The time step is set as 0.5 fs. We first equilibrate the system for 0.5 ns (10 6 steps), and then apply the external driving force for 10 ns (2 × 10 7 steps) to produce heat current and compute thermal conductivity. In both periods of simulations, the NVT ensemble (the Nose-Hoover chain thermostat) is used to control the temperature at 300 K and the thermostat relaxation time is set as 0.05 ps (100 time steps). To reduce the simulation uncertainty, for each configuration, we implement three independent simulations. The instant thermal conductivity κ(t) is directly output from GPUMD and the cumulative time-averaged thermal conductivity κ ave (t) is calculated according to Eq. (1). The converged value of κ ave (t) is taken as the final value of thermal conductivity. The reported values of thermal conductivity in this work are averaged over the three independent runs with a standard deviation. More details about the HNEMD simulation can be found in Supplementary Note 2 in Supplementary Information. The choice of driving force F e is crucial to the simulation results 69,70 . We have tested the impact of different values of F e on the thermal conductivity of graphene, which is shown in Supplementary Fig. 3 in Supplementary Information. It is seen that κ(t) converges to reasonable values when F e is about 0.05~0.1 μm −1 and behaves unexpectedly when F e > 0.1 μm −1 . Herein, we choose F e = 0.1 μm −1 to keep consistent with ref. 45,71 .
In order to eliminate the size effect in the HNEMD simulations, we carry out the convergence test with respect to the domain size. Several simulations with different domain sizes from 4.7 nm × 4.7 nm to 39.6 nm × 39.6 nm are performed. It is found that the thermal conductivity reaches a size-independent value after the size increases to about 12.3 nm × 12.3 nm (see Supplementary Note 3 in Supplementary Information). The obtained converged thermal conductivity of single-layer graphene at 300 K is 2906 ± 31 W/mK, which is very close to the value of 2903 ± 96 W/mK by EMD simulation in ref. 36 and 2847 ± 49 W/mK by HNEMD simulation in ref. 45 .
Phonon BTE. Phonon BTE is a governing equation that applies to particle transport of phonons 72 . The simulation setup of phonon BTE is shown in Fig. 1b. Periodic boundary conditions are applied at the top, bottom, left, and right boundaries of the simulation domain, which are the same as the boundary conditions applied in the HNEMD simulation. The boundary of the pore is set as a diffusely reflecting boundary 47 . To numerically solve the mode-resolved phonon BTE, the intrinsic phonon properties of materials are needed as the input. For steady-state phonon BTE, the input phonon property is the thermal conductivity accumulation distribution with respect to the phonon MFP, i.e., the cumulative MFP distribution 73 . In our previous work 42 , we reconstructed the cumulative MFP distribution of pristine graphene from NEMD simulation with Tersoff potential. Through this approach, we ensure the same intrinsic phonon properties of graphene in NEMD and BTE. A quantitative consistency of thermal conductivities predicted from NEMD and mode-resolved phonon BTE with reconstructed MFP distribution is obtained for single-layer graphene at different lengths. In this paper, the identical intrinsic phonon properties in HNEMD and BTE are also required to perform a rigorous comparison between HNEMD and BTE. In MD simulation, the intrinsic phonon properties are naturally determined by the interatomic potentials 74,75 . Since we use the same interatomic potential as the previous work 42 , we believe that the reconstructed cumulative MFP distribution can also be applicable in this work. Therefore, we still adopt this as the input of phonon BTE. The reconstructed MFP distribution is shown in Fig. 7b in ref. 42 . With the reconstructed MFP distribution, BTE predicts the thermal conductivity of pristine graphene as 2886 W/mK. This is very close to the result of the HNEMD simulation of 2906 ± 31 W/mK. The finite volume method and the discrete ordinate method are used to numerically solve the mode-resolved phonon BTE with in-house codes 76,77 . The spatial domain is discretized into unstructured grids. The angular domain at each spatial point is discretized into 128 non-overlapping control angles.

Data availability
Data supporting the findings of this study are available from the corresponding author on request.

Code availability
Computer code written and used in the analysis is available from the corresponding author as per requested.