How Hydrodynamic Phonon Transport Determines the Convergence of Thermal Conductivity in Two-Dimensional Materials

The phonon Boltzmann transport equation combined with first-principles calculation has achieved great success in exploring the lattice thermal conductivity (κ) of various materials. However, the convergence of the predicted κ is a critical issue, leading to quite scattered results recorded in the literature, even for the same material. In this paper, we explore the origin for the convergence of thermal conductivity in two-dimensional (2D) materials. Two kinds of typical 2D materials, graphene and silicene, are studied, and the bulk silicon is also compared as a control system for a three-dimensional material. The effect of the cutoff radius (rc) in the third-order interatomic force constants on κ is studied for these three materials. It is found that that κ of these three materials exhibits diverse convergence behaviors with respect to rc, which coincides very well with the strength of hydrodynamic phonon transport. By further analyzing the phonon lifetime and scattering rates, we reveal that the dominance of the normal scattering process gives rise to the hydrodynamic phonon transport in both graphene and silicene, which results in long-range interaction and a large lifetime of low-frequency flexural acoustic phonons, while the same phenomenon is absent in bulk silicon. Our study highlights the importance of long-range interaction associated with hydrodynamic phonon transport in determining the thermal conductivity of 2D materials.

The phonon Boltzmann transport equation (BTE) is a powerful tool to predict the material's lattice thermal conductivity (κ) and understand the underlying phonon transport mechanism [31][32][33][34][35][36][37][38]. Its prediction accuracy depends on a number of input parameters [39][40][41][42][43]. Among them, the cutoff radius (r c ) for the nearest neighbors' (NNs) interaction is critical in predicting the κ of 2D materials. For instance, Qin et al. [40] found that the predicted κ values for graphene by BTE are quite different with different r c , and a large r c is required to ensure the convergence of κ. Similarly, Sun et al. [44] found in four kinds of stable 2D honeycomb structures that the predicted value of κ varies greatly with different r c , which can differ by a factor of 2 if an insufficiently large r c is used, and a large r c Nanomaterials 2022, 12, 2854 2 of 13 up to the 14 th NN should be used to reach convergence. In contrast, the predicted κ of three-dimensional (3D) materials is less sensitive to the choice of r c [45,46].
Possible origins for the convergence issue for κ of 2D materials have been identified. One is that the existence of long-range interatomic forces [40,47], such as resonant bonds, which were first found in rock-salt-like crystal structures, and the decoupling of the π-bond. Another is that the truncation errors in existing calculations can lead to the divergence of κ. Lindsay et al. [48] reported that a small supercell or a finite number of nearest neighbor atomic layers may break the principle of the translational invariance (TI) condition of the crystal. Similarly, Taheri et al. [35] found that the hybridization of flexural acoustic (ZA) modes and transverse/longitudinal acoustic (TA/LA) modes in 2D materials, as well as the limitation of the supercell size could break the rotational symmetry in second-order interatomic force constants (IFCs). However, the supercell used in previous works [40,44,49] is not sufficiently large, so that the r c for the NNs' interaction exceeds half of the supercell length, which can greatly affect the prediction accuracy. Moreover, the influence of the NNs in the third-order IFCs, as well as the different sensitivity to the NNs for different dimensional systems are still unknow.
In this work, we investigate the relationship between κ and r c of the third-order IFCs in different materials, including graphene, silicene, and silicon. By calculating the value of κ with different r c , we find that the κ of graphene is the most difficult one to converge as r c increases, while silicon's κ is not sensitive to the NNs and can converge within a small r c . Through the study of phonon-phonon scattering process, we discovered that the dominant normal scattering enables long-range interatomic force in hydrodynamic 2D materials, such as graphene. By further analysis of each phonon mode, we found that the underlying physical mechanism lies in the strong anharmonic ZA phonons, while TA and LA are insensitive to r c . Our work provides physical insights into the deeper understanding of heat transport and hydrodynamic phonon transport in different systems.

Computational Method
All first-principles calculations were performed based on density functional theory (DFT), as implemented in the Vienna ab initio Simulation Package (VASP) [50][51][52] with the projected augmented wave method. The Perdew-Burke-Ernzerhof (PBE) of generalized gradient approximation (GGA) was used as the exchange-correlation functional in our study, since it has already been shown to have an accurate description of a smoother potential, correct behavior under uniform scaling, and a linear response of the uniform electron gas [53]. A vacuum space of 20 Å was imposed along the c-axis only for 2D materials to exclude the interaction between adjacent imaging layers. In all three structures, the unit cell was optimized with a cutoff energy of 600 eV until the energy and the Hellmann-Feynman force converged to 10 −8 eV and 10 −4 eVÅ −1 , respectively. Besides, the k-mesh grid used to sample the Brillouin zone (BZ) was set as 21×21×1 for 2D graphene and silicene and 15×15×15 for 3D silicon, respectively. The second-order IFC and the third-order IFCs were recorded by the finite displacement method as implemented in the PHONOPY code [54] and thirdorder_vasp.py code [55], respectively. We built a 5×5×1 supercell for graphene and silicene and a 3×3×3 supercell for silicon to obtain the second-order IFCs. In order to ensure a large NN in graphene, when calculating the third-order IFCs, we used a 7×7×1 supercell so that a large r c up to the 15 th NN can be correctly modeled, while for another two materials, the 5×5×1 supercell and 4×4×4 supercell were used for silicene and silicon, respectively. We used a finite displacement of 0.01 Å in the PHONOPY package, which is sufficient to ensure a convergent dynamical matrix.
Previous studies revealed the important role of the four-phonon scattering process in determining lattice thermal conductivity [20,21,56]. However, due to the complexity in modeling four-phonon scattering, only a short cutoff distance up to the second nearest neighbors can be considered in the calculation of fourth-order interatomic force constants [57][58][59][60]. Since our objective was to study the effect of a long-range cutoff distance in computing thermal conductivity, we limited our study to only the three-phonon scattering process.
By iteratively solving the BTE in the framework of three-phonon scattering, the lattice thermal conductivity tensor κ αβ can be obtained as implemented in the ShengBTE package [55]: where λ, , ω, k B , T, Ω, f 0 , and υ denote the phonon mode index, reduced Planck constant, phonon frequency, Boltzmann constant, temperature, system volume, equilibrium Bose-Einstein distribution, and phonon group velocity, respectively. Here, α and β are the Cartesian coordinate directions and the term F β λ denotes the linear coefficient of the nonequilibrium phonon distribution function given by where τ 0 λ is the phonon lifetime obtained under single-mode relaxation time approximation (RTA). The term ∆ β λ is obtained by the fully iterative solution of the BTE, which is a correction term to reflect the phonon distribution deviation from the RTA scheme. In the calculations, one can start with the RTA solution and obtain F β λ by iteratively solving the BTE. Based on Matthiessen's rule, the phonon lifetime τ 0 λ can be expressed as where N is the total number of q points and Γ with the superscripts (±) represents the scattering rates for three-phonon processes given by: The scattering matrix V ± λλ 1 λ 2 is given by where Φ αβγ ijk is the third-order IFCs, i, j, and k are the atomic indices, m denotes the atomic mass, and e iα λ is the α th component of the phonon eigenvector on atom i. More details about the BTE calculations can be found in our previous works [19,38].

Results and Discussion
The atomic structures of graphene and silicene are shown in Figure 1a,b, respectively. Although silicene has a similar 2D honeycomb structure as graphene, there exists a buckling height ∆ = 0.45 Å in silicene. The calculated lattice constants for graphene and silicene are Nanomaterials 2022, 12, 2854 4 of 13 2.47 Å and 3.86 Å, respectively, in good agreement with the literature results [41,43,61]. In comparison, Figure 1c shows the atomic structure of 3D bulk silicon, with a lattice constant of 5.47 Å, which also agrees well with the previous study [55]. The phonon dispersion ( Figure S1 in the Supporting Information) reveals that the quadratic dispersion of the ZA mode is only observed in 2D graphene and silicene, while all acoustic branches are linear in 3D bulk silicon.

Results and Discussion
The atomic structures of graphene and silicene are shown in Figure 1a,b, respectively. Although silicene has a similar 2D honeycomb structure as graphene, there exists a buckling height Δ = 0.45 Å in silicene. The calculated lattice constants for graphene and silicene are 2.47 Å and 3.86 Å, respectively, in good agreement with the literature results [41,43,61]. In comparison, Figure 1c shows the atomic structure of 3D bulk silicon, with a lattice constant of 5.47 Å, which also agrees well with the previous study [55]. The phonon dispersion ( Figure S1 in the Supporting Information) reveals that the quadratic dispersion of the ZA mode is only observed in 2D graphene and silicene, while all acoustic branches are linear in 3D bulk silicon. We first verified the calculated of graphene versus the Q-grid size ( × × 1). Since a 7 × 7 × 1 supercell was used for graphene in our calculations, the largest should be less than half of the length of the supercell due to the periodic boundary condition. This means the upper limit of is 8.645 Å, which corresponds to the 15 th NN ( = 8.61 Å). The detailed relation between and the n th NN can be found in Figure S2 in the Supporting Information. Figure 2a shows the convergence for the calculated of graphene at 300 K. When the cutoff radius is small (NNs = 10 th , 11 th , and 12 th ), the value of diverges with the increasing Q-grid. A similar divergent κ was also reported by Qin and Hu [40] when the is small. However, we noticed our calculation results are quantitatively different from the results of Qin and Hu [40], in which a converged κ for NNs = 10 th was reported by Qin and Hu's work (empty hexagon in Figure 2a). A 5 × 5 × 1 supercell was used in their study, which sets the upper limit of around 6.17 Å. However, = 6.35 Å (NNs = 9 th ) and = 6.80 Å (NNs = 10 th ) were used in their calculations. This discrepancy between our result and the previous study highlights the importance of using a sufficiently large supercell when increasing the cutoff radius. We first verified the calculated κ of graphene versus the Q-grid size (N×N×1). Since a 7×7×1 supercell was used for graphene in our calculations, the largest r c should be less than half of the length of the supercell due to the periodic boundary condition. This means the upper limit of r c is 8.645 Å, which corresponds to the 15 th NN (r c = 8.61 Å). The detailed relation between r c and the n th NN can be found in Figure S2 in the Supporting Information. Figure 2a shows the convergence for the calculated κ of graphene at 300 K. When the cutoff radius is small (NNs = 10 th , 11 th , and 12 th ), the value of κ diverges with the increasing Q-grid. A similar divergent κ was also reported by Qin and Hu [40] when the r c is small. However, we noticed our calculation results are quantitatively different from the results of Qin and Hu [40], in which a converged κ for NNs = 10 th was reported by Qin and Hu's work (empty hexagon in Figure 2a). A 5×5×1 supercell was used in their study, which sets the upper limit of r c around 6.17 Å. However, r c = 6.35 Å (NNs = 9 th ) and r c = 6.80 Å (NNs = 10 th ) were used in their calculations. This discrepancy between our result and the previous study highlights the importance of using a sufficiently large supercell when increasing the cutoff radius.
When the r c further increases beyond the 12 th NN, a convergent κ with respect to the Q-grid size is observed in graphene. This is because the long-range interaction can affect the lifetime of low-frequency phonons, which correspond to acoustic phonon modes near the Γ point [40]. However, the converged κ values depend notably on the specific choice of r c . For instance, the predicted κ at 300 K with the 13 th , 14 th , and 15 th NN is 5780 Wm −1 K −1 , 7300 Wm −1 K −1 , and 4867 Wm −1 K −1 , respectively, exhibiting a non-monotonic dependence on r c . A further increase of r c is not possible unless an even larger supercell than 7×7×1 is used, which is beyond our computational capacity. As shown in Figure 2a, the predicted κ value for the largest r c considered in this study (NNs = 15 th ) agrees reasonably well with the previous theoretical and experimental results [4,[62][63][64]. This complex convergence behavior in graphene is intriguing and also demonstrates that the cutoff radius is very important in determining the κ of graphene.  As another 2D material, the non-monotonic dependence of κ on r c is also observed in silicene, as shown in Figure 2b. However, when the r c is beyond the 6 th NN, the predicted κ value at 300 K in silicene reaches convergence much faster compared to that in graphene. In silicene, the converged κ at the large Q-grid limit for the 6 th , 7 th , and 8 th NN is 22.77 Wm −1 K −1 , 22.85 Wm −1 K −1 , and 21.42 Wm −1 K −1 , respectively. These values agree quite well with the results recorded in the literature [65][66][67][68]. In contrast, when switching to a 3D material, this significant convergence issue with respect to r c observed in the 2D material vanishes in the bulk silicon (Figure 2c). In the large Q-grid limit, the predicted κ value at 300 K in silicon for the 3 rd and 4 th NN is 150 Wm −1 K −1 and 154 Wm −1 K −1 , respectively, which is also in accordance with the literature results [39,[69][70][71][72][73]. By comparing these three materials, we found that the κ of the 3D material converges more easily than the 2D material, while the κ of graphene is more difficult to converge in the 2D material. This difference in the convergence of κ in the above-mentioned three materials aroused our interest in the intrinsic phonon transport form.
The different convergence behaviors of κ between the 2D and 3D materials remind us of the unique ZA branches in 2D materials. As a special feature of 2D materials, the flexural ZA phonons not only lead to unique phonon transport regimes, but also facilitate the momentum-conserving normal phonon-phonon scattering (N-scattering) [30]. When the N-scattering dominates over other phonon scatterings, such as Umklapp scattering (U-scattering), hydrodynamic phonon transport will take place [27,29].
In order to quantitatively compare the degree of hydrodynamic phonon transport among these three materials, we computed the ratio between two thermal conductivity predictions, i.e., the iterative solution (κ I ) and relaxation time approximation (κ RTA ). The relaxation time approximation (RTA) incorrectly treats the N-scattering as a resistive process, the same as U-scattering. As a result, RTA overestimates the phonon-phonon scattering rates and, thus, leads to an underestimation of κ compared to the iterative solution, especially when the N-scattering is dominant [26,28]. Therefore, the difference between κ I and κ RTA can serve as an indicator of hydrodynamic phonon transport. Figure 3 shows the ratio κ I /κ RTA versus temperature for three materials. For both graphene and silicene, this ratio is larger than unity, indicating the emergence of hydrodynamic phonon transport in these 2D materials. With decreasing temperature, this ratio gradually increases in both graphene and silicene, due to the suppressed U-scattering at a low temperature. Moreover, this ratio in graphene is larger than that in silicene at a given temperature, and the decreasing trend of this ratio with decreasing temperature is much more obvious in graphene. These results suggest a stronger hydrodynamic behavior occurring in graphene. In contrast, this ratio in bulk silicon is close to unity and is insensitive to temperature, confirming that there is no hydrodynamic phonon transport in 3D silicon.
In order to understand the origin for the different degrees of hydrodynamic phonon transport, we further calculated the phonon lifetime from both iterative and RTA solutions for these three materials. In the 2D system, there exists notable differences in the phonon lifetime between the iterative and RTA solutions for both graphene (Figure 4a) and silicene (Figure 4b), especially for the low-frequency phonons near the Γ point. These low-frequency phonons actually correspond to the ZA branch. Moreover, graphene has a mirror reflection symmetry due to the atomically flat surface, while silicene has a buckled surface. The mirror reflection symmetry in graphene imposes a symmetry selection rule on ZA phonons, causing the three-phonon processes involving an odd number of flexural phonons to be forbidden [63,74,75]. This unique selection rule in graphene gives rise to a surprisingly large phonon lifetime for ZA phonons near the Γ point, which is larger than that in silicene by more than three-orders of magnitude. In contrast, there exists a negligible difference in the phonon lifetime between the iterative and RTA solutions for the bulk silicon (Figure 4c). These results reveal that the ZA phonon is responsible for the hydrodynamic phonon transport observed in 2D graphene and silicene. In order to understand the origin for the different degrees of hydrodynamic phonon transport, we further calculated the phonon lifetime from both iterative and RTA solutions for these three materials. In the 2D system, there exists notable differences in the phonon lifetime between the iterative and RTA solutions for both graphene (Figure 4a) and silicene (Figure 4b), especially for the low-frequency phonons near the Γ point. These low-frequency phonons actually correspond to the ZA branch. Moreover, graphene has a mirror reflection symmetry due to the atomically flat surface, while silicene has a buckled surface. The mirror reflection symmetry in graphene imposes a symmetry selection rule on ZA phonons, causing the three-phonon processes involving an odd number of flexural phonons to be forbidden [63,74,75]. This unique selection rule in graphene gives rise to a surprisingly large phonon lifetime for ZA phonons near the Γ point, which is larger than that in silicene by more than three-orders of magnitude. In contrast, there exists a negligible difference in the phonon lifetime between the iterative and RTA solutions for the bulk silicon (Figure 4c). These results reveal that the ZA phonon is responsible for the hydrodynamic phonon transport observed in 2D graphene and silicene. We demonstrated the convergence issue of κ and hydrodynamic phonon transport in 2D materials, while the same phenomenon is absent in 3D materials. In order to understand the physical origin for this difference, we finally computed the frequency-resolved phonon scattering rate from Fermi's golden rule for both N-scattering (Γ N (ω)) and U-scattering (Γ U (ω)) in these three materials at 300 K. To further compare the ensemble-averaged scattering strength, the averaged scattering rate is computed as [27][28][29] where the subscripts denote the N-or U-scattering and C(ω) is the mode specific heat given by Figure 5a,b show the computed Γ N (ω) and Γ U (ω) with different NNs in graphene, respectively. For a given NN, the Γ N (ω) is much larger than Γ U (ω) by two-orders of magnitude in the low-frequency region, which makes the most contribution to the thermal conductivity. The ensemble-averaged scattering rate for the acoustic phonons is Γ N = 2.843 ps −1 and Γ U = 0.015 ps −1 . In other words, Γ N accounts for 99.5% of the total scatterings, demonstrating the dominance of N-scattering over U-scattering. This gives rise to the hydrodynamic phonon transport in graphene, as witnessed in Figure 4. In the hydrodynamic transport regime, phonons experience little scattering inside the material and, thus, can propagate collectively, which leads to very long phonon mean free path (MFP) and fast transport of thermal energy [27]. Such a long MFP and collective motion of phonons result in a long-range interaction between atoms in graphene. This effect can be revealed by the obvious dependence of Γ N (ω) and Γ U (ω) on the cutoff radius (NNs). Therefore, the hydrodynamic-transport-induced long-range interaction is responsible for the significant convergence issue of κ for graphene observed in Figure 1a.   transport of thermal energy [27]. Such a long MFP and collective motion of phonons result in a long-range interaction between atoms in graphene. This effect can be revealed by the obvious dependence of Γ ( ) and Γ ( ) on the cutoff radius (NNs). Therefore, the hydrodynamic-transport-induced long-range interaction is responsible for the significant convergence issue of κ for graphene observed in Figure 1a. Similar to graphene, Figure 5d,e reveal that the Γ N (ω) in silicene is also larger than Γ U (ω) in the low-frequency region, and ensemble-averaged scattering rate Γ N accounts for 76.9% of the total scatterings, which is larger than Γ U . However, the difference between Γ N and Γ U in silicene is much smaller than that in graphene. Furthermore, the scattering rate spectrum in silicene is less sensitive to the choice of NNs compared to that in graphene. These results suggest a weaker degree of hydrodynamic phonon transport in silicene compared to graphene. This in turn leads to a faster convergence of κ at a smaller NN ( Figure 2b) and a smaller difference in κ between the two BTE solutions (Figure 4) in silicene, compared to the case of graphene. In contrast, when it shifts to a 3D material, Figure 5g,h show no difference in the phonon scattering rate spectrum between different NNs in the bulk silicon, and the averaged scattering rate in Figure 5i confirms that the U-scattering is the dominant scattering process. This results in the vanishing of hydrodynamic phonon transport ( Figure 4) and, consequently, no convergence issue of κ for the bulk silicon, even at very small NNs (Figure 2c).

Conclusions
In summary, we calculated the room temperature thermal conductivity κ of graphene, silicene, and bulk silicon within different cutoff radius r c . In contrast to the fast convergence of κ at very small r c in the 3D bulk silicon, the κ of the 2D materials exhibits a persistent dependence on r c in both graphene and silicene, which is more pronounced in graphene even up to the 15 th nearest neighbor. By comparing the κ values predicted from the iterative solution and relaxation time approximation, the hydrodynamic phonon transport was confirmed in both graphene and silicene, while it is absent in the bulk silicon. More interestingly, the sensitivity of κ with r c in these three materials correlates very well with the strength of the hydrodynamic transport behavior. The calculations of phonon lifetime revealed that the presence of a unique ZA phonon in 2D materials is a key factor for the hydrodynamic phonon transport observed in graphene and silicene. By further analyzing the scattering rates, we demonstrated that the dominance of the N-scattering process causes the hydrodynamic transport, so that phonons can propagate collectively, which results in a long-range interaction in a 2D material such as graphene compared to the 3D material. Therefore, the long-range interaction induced by the strong hydrodynamic phonon transport is responsible for the significant convergence behavior of κ observed in graphene. This work uncovers the underlying physical mechanism for the notable convergence behavior of κ in 2D systems.

Data Availability Statement:
The data presented in this study are available from the corresponding author upon reasonable request.