Intervalley scattering induced significant reduction in lattice thermal conductivities for phosphorene

Abstract

The thermal transport properties of buckled phosphorene (β-P) and antimonene (β-Sb) are investigated using first-principles methods. The large acoustic–optical phonon gaps of 3.8 THz and 2.2 THz enable the four-phonon interaction to play an important role in phonon scattering for both β-P and β-Sb. Considering the electron–phonon coupling, the lattice thermal conductivity can further undergo 84% decrease to 4.9 W mK⁻¹ for p-type β-P at n = 5 × 10¹³ cm⁻². By quantitatively describing the scattering probability of electrons in different paths combined with electron–phonon coupling matrix element analysis, it is found that multi-valley features of electronic band structure and strong electron–phonon coupling strength make electrons have strong intervalley scattering behavior in β-P. The former plays an important role in the energy conservation condition of the scattering process, and the latter determines the selection rule. Our work elucidates the contribution of higher-order phonon interactions as well as electron–phonon coupling effects to lattice thermal conductivity, and provides a new idea for finding materials with low lattice thermal conductivity induced by intervalley scattering.

---

New concepts

The interaction between phonons has long been considered as the dominant factor in phonon scattering in materials. However, we find that the intervalley scattering effect caused by the electron–phonon coupling dominates the phonon scattering in phosphorene with multi-valley features, and the interaction of four phonons in the material cannot be ignored due to the large acoustic–optical phonon gap. Our work provides a new idea for finding materials with low lattice thermal conductivity induced by intervalley scattering.

1 Introduction

When the carrier concentration is not very high, the intrinsic phonon–phonon interaction usually dominates the phonon scattering.^1,2 Electron–phonon coupling has long been considered to limit the electrical transport performance and has a trivial effect on phonon thermal transport in lowly doped semiconductors. The lattice thermal conductivities κ_L calculated by the first-principles method considering phonon–phonon interaction can be in good agreement with the experimental results in a wide variety of materials with low carrier concentration.^3–6 However, with the improvement of doping technology, the carrier concentration can reach 10²¹ cm⁻³ for 3D materials and 10¹⁴ cm⁻² for 2D materials,^7,8 motivating researchers to reconsider the contribution of electron–phonon coupling to phonon scattering.^9–11

For instance, Liao et al. reported a dramatic reduction of lattice thermal conductivity by as much as 45% for bulk silicon considering the electron–phonon coupling effect with 10²¹ cm⁻³ carrier concentration at room temperature.⁹ Li et al. found that the indirect interaction between the ZA phonons and electrons mediated by the TA and LA phonons in graphene can result in 21% reduction of lattice thermal conductivity at room temperature.¹² Liu et al. explained the abnormal suppression effect of electron–phonon interaction on lattice thermal conductivity in p-type PtSSe by analyzing the electron–phonon coupling strength and electronic band characteristics.¹³

The dispersion of electronic structure plays an important role in the energy conservation conditions of the electron–phonon coupling process and then further affects the phonon scattering rate. In general, the higher the density of states near the Fermi level, the more scattering channels of electrons, and the stronger the electron–phonon interaction. Group-V elemental materials (P, As, Sb, Bi) have a stable buckled honeycomb structure and have been successfully synthesized via molecular beam epitaxy (MBE) or pulsed laser deposition (PLD).^14–17 They have diverse band structures and prove to own many intriguing and useful electrical transport properties.^18–22 Previous calculations show that β-P has ultra-low hole mobility (0.9 cm² V⁻¹ s⁻¹) due to the strong electron–phonon coupling effect induced by its extremely flat valence band structure with multi-valleys.²³ We expect that this strong coupling effect can in turn efficiently increase the phonon scattering rate and limit the phonon transport properties. Here, the effect of electron–phonon coupling on the lattice thermal conductivity of β-P is systematically investigated from first principles calculations and compared with the results of β-Sb with the same structure. The four-phonon scattering effect is included due to the large acoustic–optical (a–o) gap in both β-P and β-Sb, which is considered as the inducement of strong four-phonon scattering.^24,25 The results show that four-phonon scattering and the intervalley scattering effect caused by electron–phonon coupling play important roles in phonon scattering and can effectively reduce the lattice thermal conductivity.

2 Results and discussion

The configurations of the monolayer β-P/Sb are shown in Fig. S1 (ESI†). Different from the planar geometry of graphene, the monolayer β-P/Sb has a buckled structure with two sublayers and thus the horizontal mirror symmetry is broken. The broken mirror symmetry can lead to different phonon–phonon and electron–phonon coupling strengths from graphene. Actually, studies have shown that symmetry breaking can significantly reduce the phonon degeneracy and enhance the avoid-crossing behavior of phonon branches, which results in an increase in the phonon–phonon scattering phase space to limit lattice heat transport.^26,27 The optimized lattice constants are 3.28 Å and 4.12 Å for β-P and β-Sb, respectively, which are in good agreement with the previous work.^28,29 The calculated band structures are shown in Fig. 1(a and b) and the corresponding energy of the center of the band gap is set to zero. The indirect band gaps of β-P and β-Sb are 1.94 eV and 0.99 eV, respectively. The conduction band minimum (CBM) of β-P and β-Sb locates at the C₁ point close to the middle of the Γ–M direction. It is noted that the band structure of β-P near the CBM is particularly flat, which will significantly enhance the electron–phonon coupling. The valence band maximum (VBM) of β-P locates at the V₁ point along the Γ–K direction, while for β-Sb the VBM locates at the Γ point. There exists high band degeneracy in the Brillouin Zone for both systems, which may induce intervalley scattering channels of electrons in the electron–phonon coupling process. The band degeneracy near the CBM is N_c = 6 for both β-P and β-Sb. The valence band of β-P possesses local peaks V₂ (N_v = 6) and V₃ (N_v = 2) below the global VBM (N_v = 6) that will play an important role in the lattice thermal conductivity of the system under p-type doping. The energy difference between V₂ and V₁ (V₃ and V₁) is 49 meV (162 meV). The insets of Fig. 1(a and b) show the 3D band structure of the valence band for β-P and β-Sb, respectively. The valence band of β-P possesses multiple valleys and peaks which are shallow in energy and connected with each other. In contrast, only a sharp and deep peak locates at the Γ point for β-Sb and the holes can only experience intravalley scattering in the electron–phonon coupling process.

Fig. 1 The electronic band structure of (a) β-P and (b) β-Sb. C₁ and V₁ are the location of the CBM and VBM. V₂ and V₃ denote local band-maxima in the valence band of β-P. The inset is the 3D band structure of the valence band. Phonon dispersion for (c) β-P and (d) β-Sb with the color representing the average phonon occupancy at 300 K.

The phonon dispersion curves of β-P and β-Sb are shown in Fig. 1(c and d) with the color of the curve indicating the average phonon occupancy at room temperature calculated using the Bose-Einstein distribution 1/(e^hω/k_BT − 1). Different from the planar graphene, the buckled structure of β-P and β-Sb results in the coupling of in-plane and flexural phonons, leading to a considerable a–o gap of 3.8 THz and 2.2 THz, respectively. The large a–o gap makes the excited phonons mainly concentrated on the acoustic branch and the phonon occupation number in the optical branch is mainly less than unity. In general, the large a–o gap is associated with significant four-phonon scattering due to the fact that the total energy of three acoustic phonons is large enough to cross the a–o gap recombining into an optical phonon.

The calculated T-dependent lattice thermal conductivity κ_L is shown in Fig. 2(a and b). The calculated κ_L of 114.1 W mK⁻¹ and 10.7 W mK⁻¹ at 300 K for β-P and β-Sb under the perspective of three-phonon scattering is close to the previous calculations.^30,31 We find that the four-phonon and electron–phonon interaction have strong effects on the lattice heat transport properties of β-P and β-Sb. Specifically, at 300 K, the κ_L is 30.1 W mK⁻¹ (74% decrease) for β-P when including the four-phonon scattering in the pristine case and further reduces to 4.9 W mK⁻¹ (84% decrease) and 14.7 W mK⁻¹ (51% decrease) by p-type and n-type doping, respectively. For β-Sb, the κ_L is 2.0 W mK⁻¹ (81% decrease), 1.6 W mK⁻¹ (20% decrease) and 1.3 W mK⁻¹ (35% decrease) when considering four-phonon scattering for the pristine, p-type and n-type system, respectively. This indicates that four-phonon interaction dominates the phonon scattering in both β-P and β-Sb while the effect of electron–phonon coupling on phonon scattering is more obvious in β-P. Moreover, it can be found that the κ_L with four-phonon scattering exhibits stronger temperature dependence while the case is the opposite when the electron–phonon coupling is taken into account. The κ_L is nearly inversely proportional (∝T⁻¹) to temperature when only considering three-phonon scatterings, which is due to the dominance of Umklapp processes among the total scatterings. At high temperature, more phonons are excited, and the four-phonon interaction is more obvious, making κ_L more dependent on temperature. On the contrary, the phonon scattering from electron–phonon coupling is almost T-independent. Therefore the impact of electron–phonon coupling on κ_L is more obvious at lower temperature. As seen in Fig. 2(c, d), considering phonon–phonon scatterings, the κ_L of β-P and β-Sb is mainly dominated by acoustic phonons, and the optical phonons have less contribution to thermal transport. Therefore, the contribution of electron–phonon coupling effect on phonon scattering only needs to consider the participation of acoustic phonons. The insets of Fig. 2(c and d) show the contribution of different acoustic phonon modes to κ_L. LA modes play a very important role in thermal transport when the effect of electron–phonon coupling is not considered. This is due to the large a–o gap in β-P and β-Sb, which suppresses the recombination of LA phonons into optical phonons. In order to clearly describe the phonon recombination process, the phase space of three phonon interaction for recombination process P3⁺ is calculated and shown in Fig. S2 (ESI†). For β-P, in the phonon frequency of 0–4 THz, the P3⁺ of the LA phonons is mainly lower than other acoustic branches. When the frequency is further increased, some of the LA phonons can be recombined into optical phonons, leading to gradual increases of P3⁺. For β-Sb, since the maximum frequency of LA phonons is smaller than the a–o gap, the process of LA phonon recombination into optical phonons is completely inhibited, and the P3⁺ of the LA phonons remains small in a wide frequency range. However, the inclusion of the electron–phonon coupling makes TA phonons become the main source of thermal transport.

Fig. 2 The variation of lattice thermal conductivity κ_L with temperature considering different sources of phonon scattering for (a) β-P and (b) β-Sb. The calculated result considering three-phonon interaction for β-P from reference is also shown with green points (ref. 30). Accumulation functions of κ_L considering phonon–phonon interaction at 300 K for (c) β-P and (d) β-Sb. The insets show the contribution of ZA, TA, LA and optical phonon modes to κ_L.

The κ_L for β-P and β-Sb at 300 K, varying with carrier concentration, is shown in Fig. 3(a and b). It is found that the variation of κ_L with carrier concentration generally presents a trend of rapid decrease first and then a slow increase. As a result, the κ_L reaches its minimum value near n = 2.5 × 10¹⁴ cm⁻² (n = 4.5 × 10¹⁴ cm⁻²) and n = 2.5 × 10¹⁴ cm⁻² (n = 3.5 × 10¹⁴ cm⁻²) for p-type (n-type) β-P and p-type (n-type) β-Sb, respectively. The electrons participating in the electron–phonon coupling are mainly concentrated near the Fermi level. Hence, when more electrons occupy the state near the Fermi level, the phonons are more strongly scattered from electron–phonon coupling. Carrier concentration is a function of the Fermi level E_F, defined as^32,33

where Ω is the unit cell area, f₀ is the Fermi-Dirac distribution, ε is the energy of carriers and D(ε) is the density of states (DOS). Since a higher carrier concentration is required when the electron–phonon coupling effect has a significant impact on phonon transport, making the Fermi level outside the band gap, we can calculate the corresponding DOS and n at a given Fermi level. Fig. 3(c and d) show the electron and hole DOS of β-P and β-Sb with respect to the carrier concentration n. The carrier concentration n as a function of Fermi energy E_F is shown in Fig. S3 (ESI†). It can be seen that the n corresponding to the maximum value of the DOS is located near the minimum of κ_L. The electron DOS of β-P oscillates around 0.5 eV⁻¹ when the n rises to a certain value, making the κ_L for n-type β-P exhibit saturation effect.

Fig. 3 The lattice thermal conductivity versus the carrier concentration n for doped (a) β-P and (b) β-Sb at 300 K. Electron density of states of (c) β-P and (d) β-Sb with respect to the carrier concentration n.

In addition to the lattice thermal conductivity κ_L, carriers can also contribute to heat transport, denoted as electronic thermal conductivity κ_e. The κ_e exhibits a similar trend with electronic conductivity σ according to the Wiedemann–Franz law,³⁴i.e. κ_e = LσT, where L is the Lorenz number. Hence, increasing σ leads to the increase of κ_e. A stronger electronic conductivity often indicates that the carriers are less scattered by phonons, and the electron–phonon coupling is weaker, which will have less influence on the lattice thermal conductivity. According to the Boltzmann transport theory, Fig. S4(a and b) (ESI†) show the κ_eversus carrier concentration of β-P and β-Sb at 300 K. The κ_e remains small over the range of considered concentration for β-P due to the strong electron–phonon coupling induced by intervalley scattering. Compared with κ_L, the contribution of κ_e to the total thermal conductivity is almost negligible in β-P. For β-Sb, the κ_e is comparable to κ_L for the n-type system. Due to the large dispersion of the valence band in β-Sb, κ_e dominates the thermal transport for the p-type system. Fig S4(c and d) (ESI†) show the total thermal conductivities κ_total for β-P and β-Sb. It can be seen that the profile of κ_total for β-P has hardly changed compared with κ_L. However, for β-Sb, the value and trend of κ_total have changed significantly.

For more insights, Fig. 4 shows the phonon scattering rates of β-P and β-Sb contributed from phonon–phonon interaction and electron–phonon coupling at 300 K with n = 5 × 10¹³ cm⁻². The value of the phonon–phonon scattering rate is small in the low frequency region and increases slowly with the increase of the frequency. This is because the phonon emission process is limited in the low frequency region. The phonon scattering rate coming from the electron–phonon coupling does not show obvious frequency dependence. The four-phonon scattering rate is comparable to the three-phonon scattering rate in β-P and even exceeds the three-phonon scattering rate in the high frequency region in β-Sb. This reveals the importance of four-phonon scattering in materials with large a–o gap. The LA phonons possess higher group velocity due to large phonon dispersion and smaller scattering rate under phonon–phonon interaction, making it contribute more to κ_L. It is clear that the electron–phonon coupling plays a dominant role in phonon scattering of β-P and covers a broad frequency range. In particular, the scattering rate of the ZA phonon reaches ∼100 ps⁻¹, three orders of magnitude higher than the value under phonon–phonon interaction. Compared with TA phonons, the scattering rate of ZA and LA phonons from electron–phonon coupling is much higher than the three-phonon and four-phonon scattering rate, resulting in TA phonons dominating thermal transport in Fig. 2(c). For β-Sb, the LA phonons in the frequency range of 0–1 THz have an electron–phonon coupling scattering rate comparable to that of phonon–phonon interactions for both p-type and n-type systems. However, a small number of ZA and TA phonons in a narrow frequency range have an electron–phonon coupling scattering rate with a significant value.

Fig. 4 The three-phonon (3ph), four-phonon (4ph), and electron–phonon (ep) scattering rates at n = 5 × 10¹³ cm⁻² of (a, d) ZA, (b, e) TA and (c, f) LA modes at T = 300 K for β-P and β-Sb.

To explain the large difference in the electron–phonon coupling scattering rate, we define a physical quantity P reflecting the scattering probabilities of electrons in the electron–phonon coupling process, written as,

where |g|² reflects the selection rule of electron–phonon coupling and determines the intensity of electron–phonon coupling, and weight includes the energy conservation and momentum conservation conditions in the process of electron–phonon coupling determining the contour of electron–phonon coupling in k space. The P parameter of holes as a function of the wave vector of final state k + q involving different acoustic phonon modes is shown in Fig. 5. Here the initial state is chosen near the VBM of β-P and β-Sb represented by an orange star in the figure. The results of the initial state near the almost degenerate point V₂ for β-P are shown in Fig. S5 (ESI†). Due to the small energy fluctuation of the valence band in β-P, the energy conservation condition of the electron–phonon coupling process is easily satisfied, and holes have a greater scattering probability in the energy degenerate valleys. Among the three types of acoustic phonons, the ZA mode has larger P especially near the initial state and the LA mode covers a wide range in k space. When setting the initial state near V₂, scattering probabilities exhibit a similar pattern, as shown in Fig S5 (ESI†). The reason for the large P coming from coupling between electrons and ZA phonons can be explained from the perspective of energy conservation. Near the initial state, the phonons with small momentum participate in the electron–phonon coupling and the typical quadratic dispersion of the ZA mode matches well with the electronic dispersion. However, the LA mode has large energies and phonon dispersion, which should not cause large electron scattering probability. For understanding this result, Fig. 6 shows the electron–phonon coupling matrix element |g|. It is clear that the energy valley with larger scattering probability usually has a larger |g| value and ZA and LA modes have a wider region with larger |g| than the TA mode. For β-Sb, the valence band has a single valley and the electron–phonon coupling effect can only induce intravalley scattering. Due to the sharp and deep electronic band structure, only a circular boundary area near the Γ point has obvious scattering probability. It can be seen that LA phonons with small momentum play a major role in electron–phonon coupling mainly due to the fact that the strong linear dispersion of the LA phonon mode matches well with electronic structure. The corresponding results of the conduction band are shown in Fig. S6 (ESI†). The electrons near the CBM of β-P and β-Sb exhibit an extremely strong intervalley scattering effect and different acoustic phonons all make considerable contributions.

Fig. 5 The scattering probability P of holes at different positions in the Brillouin Zone considering electron–phonon coupling for (a)–(c) β-P and (d)–(f) β-Sb. The orange star represents the initial states of holes.

Fig. 6 The electron–phonon coupling matrix elements |g| of the holes in β-P involving (a) ZA, (b) TA and (c) LA phonons. The orange star represents the initial states of holes.

In order to further understand the influence of the electron–phonon coupling effect on the lattice thermal conductivity of these two materials, we distinguish the contribution of phonons with different momentum to κ_L considering three-phonon interactions. As shown in Fig. 7, when considering phonon–phonon interaction, phonons with small (0–0.1 Å⁻¹) and medium momentum (0.1–0.3 Å⁻¹) mainly contribute to κ_L for β-P. The former mainly participate in the intravalley scattering of electrons in the electron–phonon coupling process, and the latter mainly participate in the intervalley scattering. For β-Sb, except for phonons near 0.02 Å⁻¹ contributing much to κ_L, phonons with other momentum have a similar contribution. Considering p-type doping with n = 5 × 10¹³ cm⁻², the change of κ_L (Δκ_L/κ_L) is depicted as a green sphere in Fig. 7. It is clear that due to the special valence band structure of β-P, the electron–phonon coupling effect obviously limits the lattice thermal conductivity over a wide range of momenta. However, since only intravalley scattering is involved in electron–phonon coupling for β-Sb, the contribution of phonons with small momentum to lattice thermal conductivity is greatly affected by electron–phonon coupling.

Fig. 7 The contribution of phonons with different momentum to the lattice thermal conductivity considering three-phonon interaction for (a) β-P and (b) β-Sb. The green sphere represents the change of thermal conductivity (Δκ_L/κ_L) under p-type doping with n = 5 × 10¹³ cm⁻².

3 Conclusion

In conclusion, we investigate the lattice thermal transport properties of β-P and β-Sb through first-principles calculations. The large a–o gap makes the acoustic phonons dominate the thermal transport and the four-phonon interaction plays an important role in the phonon scattering. Considering the phonon–phonon interaction, LA phonons dominate the heat transport due to their large phonon dispersion and low scattering rate. A significant reduction of lattice thermal conductivity occurs when the electron–phonon coupling effect is included in phonon scattering, especially for p-type β-P with an 84% decrease at n = 5 × 10¹³ cm⁻². The large electronic density of states often induces a large phonon scattering rate, and the minimum of lattice thermal conductivity appears near the corresponding carrier concentration. By analyzing the electron scattering probability and electron–phonon coupling matrix, we find that a strong intervalley scattering effect exists in the multi-valley degenerate system, which can cause a sharp drop in lattice thermal conductivity. Our work provides a deep understanding of higher-order phonon interactions and electron–phonon coupling effects on lattice thermal conductivity and provides theoretical guidance for finding novel materials with low lattice thermal conductivity.

Numerical method and computational details

According to the Boltzmann transport equation, the lattice thermal conductivity tensor can be expressed as,³⁵where N_q is the number of q points in the first Brillouin zone, and C_qλ, v_qλ^α (v_qλ^β) and τ_qλ are the heat capacity, group velocity along α (β) direction and relaxation time of phonons with branch λ and wave vector q, respectively. Considering both the phonon–phonon scattering and electron–phonon scattering effect, the relaxation time can be defined by Matthiessen's law,

1/τ_qλ = 1/τ^PP_qλ + 1/τ^ep_qλ (5)

where τ^PP_qλ is the relaxation time for phonon–phonon scattering and τ^ep_qλ is the imaginary part of the phonon self-energy due to the electron–phonon scattering, written as,³⁶where f_nk and ε_nk are the Fermi-Dirac distribution function and energy for electrons with band branch n and wave-vector k, respectively. ω_qλ is the frequency of phonons. g_mnλ(k,q) is the electron–phonon coupling matrix element.

The density functional theory (DFT) calculations are carried out by the QUANTUM ESPRESSO code.³⁷ During the processes of self-consistent and non self-consistent calculations, the kinetic energy cutoff is set to be 90 Ry, and the Perdew–Burke–Ernzerhof (PBE) exchange correlation functional³⁸ and the fully relativistic norm-conserving pseudopotentials that include the spin-orbit coupling (SOC) are employed.³⁹ A fine final grid of 120 × 120 × 1 for both k and q is used for calculating τ_qλ^ep implemented in the EPW package.⁴⁰ We use BoltzTraP2 code based on the rigid band approach to calculate the electronic thermal conductivity.³³ In this code, the input of electron relaxation time is treated as both band- and momentum-dependent obtained from EPW. The phonon–phonon interaction is calculated by the ShengBTE package³⁵ with an iterative self-consistent algorithm where a 8 × 8 × 1 and 6 × 6 × 1 supercells are built for the harmonic and anharmonic interatomic force constants (IFCs), respectively. The eight and three nearest neighboring atoms are considered for third-order and forth-order interactions, respectively. We have modified the code to consider the electron–phonon coupling effect on the lattice thermal conductivity. To get converged thermal conductivities, a q-mesh of 120 × 120 × 1 and 40 × 40 × 1 is used for the calculation of the third-order and forth-order phonon relaxation time, respectively. To deal with 2D thermal conductivity, an effective thickness of 5.63 Å and 3.76 Å is adopted for β-P and β-Sb,³⁰ respectively.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work is supported by the Shanghai Municipal Natural Science Foundation under Grant No. 19ZR1402900.

References

1. P. Klemens, Thermochim. Acta, 1993, 218, 247–255.
2. J. M. Ziman and P. W. Levy, Phys. Today, 1961, 14, 64–66.
3. D. A. Broido, M. Malorny, G. Birner, N. Mingo and D. A. Stewart, Appl. Phys. Lett., 2007, 91, 231922.
4. J. Garg, N. Bonini, B. Kozinsky and N. Marzari, Phys. Rev. Lett., 2011, 106, 045901.
5. S. Li, Q. Zheng, Y. Lv, X. Liu, X. Wang, P. Y. Huang, D. G. Cahill and B. Lv, Science, 2018, 361, 579–581.
6. F. Tian, B. Song, X. Chen, N. K. Ravichandran, Y. Lv, K. Chen, S. Sullivan, J. Kim, Y. Zhou, T.-H. Liu, M. Goni, Z. Ding, J. Sun, G. A. G. U. Gamage, H. Sun, H. Ziyaee, S. Huyan, L. Deng, J. Zhou, A. J. Schmidt, S. Chen, C.-W. Chu, P. Y. Huang, D. Broido, L. Shi, G. Chen and Z. Ren, Science, 2018, 361, 582–585.
7. D. K. Efetov and P. Kim, Phys. Rev. Lett., 2010, 105, 256805.
8. C.-L. Wu, H. Yuan, Y. Li, Y. Gong, H. Y. Hwang and Y. Cui, Nano Lett., 2018, 18, 2387–2392.
9. B. Liao, B. Qiu, J. Zhou, S. Huberman, K. Esfarjani and G. Chen, Phys. Rev. Lett., 2015, 114, 115901.
10. J.-Y. Yang, G. Qin and M. Hu, Appl. Phys. Lett., 2016, 109, 242103.
11. Y. Zhang, Z. Tong, A. Pecchia, C. Yam, T. Dumitrică and T. Frauenheim, Nanoscale, 2022, 14, 13053–13058.
12. X. Yang, A. Jena, F. Meng, S. Wen, J. Ma, X. Li and W. Li, Mater. Today Phys., 2021, 18, 100315.
13. C. Liu, M. Yao, J. Yang, J. Xi and X. Ke, Mater. Today Phys., 2020, 15, 100277.
14. J. L. Zhang, S. Zhao, C. Han, Z. Wang, S. Zhong, S. Sun, R. Guo, X. Zhou, C. D. Gu, K. D. Yuan, Z. Li and W. Chen, Nano Lett., 2016, 16, 4903–4908.
15. M. Fortin-Deschênes and O. Moutanabbir, J. Phys. Chem. C, 2018, 122, 9162–9168.
16. Y. Shao, Z.-L. Liu, C. Cheng, X. Wu, H. Liu, C. Liu, J.-O. Wang, S.-Y. Zhu, Y.-Q. Wang, D.-X. Shi, K. Ibrahim, J.-T. Sun, Y.-L. Wang and H.-J. Gao, Nano Lett., 2018, 18, 2133–2139.
17. Z. Yang, Z. Wu, Y. Lyu and J. Hao, InfoMat, 2019, 1, 98–107.
18. Z. Wu and J. Hao, npj 2D Mater. Appl., 2020, 4, 4.
19. D. Kecik, V. O. Özçelik, E. Durgun and S. Ciraci, Phys. Chem. Chem. Phys., 2019, 21, 7907–7917.
20. Y. Wang, P. Huang, M. Ye, R. Quhe, Y. Pan, H. Zhang, H. Zhong, J. Shi and J. Lu, Chem. Mater., 2017, 29, 2191–2201.
21. Y. Wu, B. Hou, Y. Chen, J. Cao, H. Shao, Y. Zhang, C. Ma, H. Zhu, R. Zhang and H. Zhang, npj Comput. Mater., 2021, 7, 145.
22. Y. Wu, B. Hou, C. Ma, J. Cao, Y. Chen, Z. Lu, H. Mei, H. Shao, Y. Xu, H. Zhu, Z. Fang, R. Zhang and H. Zhang, Mater. Horiz., 2021, 8, 1253–1263.
23. L. Cheng, C. Zhang and Y. Liu, J. Am. Chem. Soc., 2019, 141, 16296–16302.
24. J. Chen, J. He, D. Pan, X. Wang, N. Yang, J. Zhu, S. A. Yang and G. Zhang, Sci. China: Phys., Mech. Astron., 2022, 65, 117002.
25. X. Yang, T. Feng, J. Li and X. Ruan, Phys. Rev. B, 2019, 100, 245203.
26. C. Liu, C. Wu, T. Song, Y. Zhao, J. Yang, P. Lu, G. Zhang and Y. Chen, ACS Appl. Energy Mater., 2022, 5, 15356–15364.
27. C. Liu, Z. Chen, C. Wu, J. Qi, M. Hao, P. Lu and Y. Chen, ACS Appl. Mater. Interfaces, 2022, 14, 46716–46725.
28. S. Balendhran, S. Walia, H. Nili, S. Sriram and M. Bhaskaran, Small, 2014, 11, 640–652.
29. D.-C. Zhang, A.-X. Zhang, S.-D. Guo and Y. F. Duan, RSC Adv., 2017, 7, 24537–24546.
30. G. Zheng, Y. Jia, S. Gao and S.-H. Ke, Phys. Rev. B, 2016, 94, 155448.
31. B. Peng, D. Zhang, H. Zhang, H. Shao, G. Ni, Y. Zhu and H. Zhu, Nanoscale, 2017, 9, 7397–7407.
32. J. Yang, H. Li, T. Wu, W. Zhang, L. Chen and J. Yang, Adv. Funct. Mater., 2008, 18, 2880–2888.
33. G. K. Madsen, J. Carrete and M. J. Verstraete, Comput. Phys. Commun., 2018, 231, 140–145.
34. X. Zheng, D. Cahill, P. Krasnochtchekov, R. Averback and J. Zhao, Acta Mater., 2007, 55, 5177–5185.
35. W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758.
36. F. Giustino, Rev. Mod. Phys., 2017, 89, 015003.
37. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
38. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
39. D. R. Hamann, Phys. Rev. B, 2013, 88, 085117.
40. J. Noffsinger, F. Giustino, B. D. Malone, C.-H. Park, S. G. Louie and M. L. Cohen, Comput. Phys. Commun., 2010, 181, 2140–2148.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3nh00090g

---