Yi-Ming
Zhao
a,
Chun
Zhang
bc,
Sunmi
Shin
*a and
Lei
Shen
*a
aDepartment of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, Singapore, 117575, Singapore. E-mail: mpeshin@nus.edu.sg; shenlei@nus.edu.sg
bDepartment of Physics, Faculty of Science, National University of Singapore, 2 Science Drive 3, Singapore, 117551, Singapore
cDepartment of Chemistry, National University of Singapore, 3 Science Drive 3, Singapore, 117543, Singapore
First published on 18th July 2023
Two-dimensional bilayer structures exhibit novel properties not existing in the monolayer ones like bilayer hexagonal boron nitride (h-BN) with exotic interfacial ferroelectrics originated by laterally sliding one layer over the other (K. Yasuda, et al., Science, 2021, 372, 1458; M. Vizner Stern, et al., Science, 2021, 372, 1462). This intriguing property provides the sliding bilayer h-BN with promising applications in slidetronics. However, the thermal transport behaviour in the sliding bilayer h-BN is still unclear. Here, we systematically investigate the influence of the lattice configurations on the thermal conductivity (κ) of sliding bilayer h-BN structures as well as its manipulation with strain and layer confinement via first-principles. Our results show that structures with boron head-to-head stacking (B–B) exhibit lower κ values than the ones with nitrogen on the top of boron stacking (B–N). The phonon spectra and weighted phase space indicate a softer layer-breathing mode (ZO′) mode and higher three-phonon scattering rate of B–B patterns, leading to a lower κ. The moderate out-of-plane compressive strain of −6% significantly decreases the κ of B–B structures by about 50% through enhanced anharmonic scattering, while the higher strain of −18% suppresses the anharmonic scattering and increases the κ instead. Finally, we study the stacking dependence of the κ in tri-layer and bulk structures.
Thermal conductivity (κ) is an important parameter for the thermal management of electric devices. Materials with higher κ are preferable for high-power devices and efficient heat dissipation,12–14 while a low κ is crucial for thermoelectric applications.15–18 The κ of monolayer h-BN reaches 751 W m−1 K−1 based on experimental measurements, which is promising for the heat dissipation applications.19 The previous computational results based on the solution of the phonon Boltzmann transport equation (BTE) reveal a trend of decreasing thermal conductivity with the increasing number of h-BN layers.20 The results from molecular dynamics (MD) simulations demonstrated that the mismatch of different layers in bulk h-BN effectively influences the out-of-plane κ.21 Gholivand and Donmezer also confirmed that the thermal conductivity of bilayer h-BN (401 W m−1 K−1, in-plane) is lower than that of single layer h-BN (579 W m−1 K−1).22 The out-of-plane acoustic (ZA) phonon modes contribute most to the thermal conductivity in h-BN. The confinement from adjacent layers effectively changes the ZA phonon modes and decreases the κ for bilayer structures. The van der Waals sliding will change the configuration of bilayer h-BN along the out-of-plane direction, intuitively influencing the thermal conductivity to a certain extent. However, the sliding effect on the thermal conductivity is not clear yet.
Moreover, strain will change phonon vibration behavior, effectively affecting the thermal transport properties.23–25 The hydrostatic pressure drives the thermal transport from anisotropic to isotropic in layered InSe via modulating the phonon group velocity and lifetime.26 Sun et al. explored the effect of interlayer interaction on both the in-plane and out-of-plane thermal conductivity of bilayer graphene and graphite.27 An out-of-plane compressive strain decreases the in-plane κ, while it increases the out-of-plane κ. Duan et al. also investigated the strain effect on thermal conductivity and revealed an enhanced hydrodynamic phonon transport in bilayer graphene with larger strains.28 In the experiment, the diamond anvil cell29 or the piston cylinder cell30,31 is usually used to apply the hydrostatic pressure and tune the out-of-plane strain on layered materials.32
In this work, we systematically investigate the sliding and strain effects on the thermal conductivity of bilayer h-BN. Our results indicate that the thermal conductivity of AC′ (B–B) (see configurations in Fig. 1) is lower than those of both AA′ (B–N) and AB (B–N) due to the softer ZO′ modes. The thermal conductivity of AC′ can be modulated by out-of-plane compressive strains which change the frequency of ZO′ modes and the anharmonic interaction. The AC′ configuration exhibits a higher electrical conductivity compared with AA′ under p-type doping. Finally, we find that the trend in sliding bilayer h-BN is general and valid in tri-layer and bulk structures as well: i.e., lower thermal conductivity in AC′ compared with AA′ and AB stacks. This work provides a fundamental understanding of the sliding effect on the thermal transport behavior of stacked h-BN, which sheds light on the design of h-BN related thermal devices with sliding and/or strain.
The convergence test is a necessary step to obtain an accurate thermal conductivity and the number of nearest neighbors is an important parameter in the calculation of third order IFC. We, thus, calculate the κ of AA′ configuration with the number of nearest neighbors varying from 8 to 20 for the convergence test as shown in Fig. 2b. As can be seen, the thermal conductivity of AA′ finally converged at 428 W m−1 K−1 when considering the 19th nearest neighbor, which is in good agreement with the published result of AA-type h-BN in Table 1. The cutoff distance for the 19th nearest neighbor is 6.97 Å larger than the distance (5.02 Å) evaluated from the variation of charge density induced by the displacement of central boron atom.45 The 0.02 Å displacement will influence the charge distribution within the range of 5.02 Å as shown in the inset Fig. 2b. The larger cutoff distance ensures the accuracy of κ calculation. Moreover, our result for monolayer h-BN is also close to those of previous works (Table 1). The converged and reasonable results indicate that our computational method and process are correct. After the convergence test, we calculate the thermal conductivity of AB, which is close to that of AA′ because the atomic stacking pattern of AB is B–N, the same as that of AA. The thermal conductivity of AC′ with the B–B pattern is the lowest one among the three configurations, 397 W m−1 K−1. The lower κ on the out-of-plane for the bulk structure indicates an anisotropic thermal transport, which is a typical characteristic for the layered structures.46,47 Because of the same B–N pattern in AB and AA′, we group AB and AA′ together in the following discussion unless otherwise stated, and only make comparisons between the configurations of AC′ (B–B) and AA′ (B–N).
To understand the underlying physics of the lower thermal conductivity of B–B type AC′ than the B–N type AA′ and AB, we investigate and compare intermediate physical quantities between AC′ and AA′, including the weighted phase space (WP3), the Grüneisen parameter (γ), and the total scattering rate in Fig. 3. The cumulative thermal conductivity of AC′ and AA′ gradually increase to the converged value when more phonon modes with higher frequencies are involved as shown in Fig. 3a. The first order derivative of the thermal conductivity (inset of Fig. 3a) reflects the contribution of phonon modes with specific frequency to the final thermal conductivity. The contributions of AC′ phonon modes are obviously less than that of AA′ in the range from 2 to 10 THz (highlighted by a pink shadow). The negative spectral thermal conductivity near 15 THz should be a numerical artifact during the iterative process as shown in Fig. S3 (ESI†). As mentioned above, the lower frequencies of ZO′ modes in AC′ are the main difference in the phonon band structures. Based on the energy and momentum conservation rule in phonon scattering behavior, the phonon band structure will directly determine the number of 3-phonon scattering processes, also called phase space. Fig. 3b shows that the phase space with frequency weights48 of AC′ is slightly higher than AA′ in the range from 2 to 10 THz, indicating a higher scattering possibility and lower thermal conductivity of AC′. Therefore, the low frequency of ZO′ modes in AC′ is an important factor for its lower thermal conductivity.
Besides the weighted phase space, other intermediate physical quantities like the γ parameter and scattering rate also contribute to the difference of the thermal transport behavior between AC′ and AA′. The γ parameter measures the anharmonicity of the lattice vibration, thereby indicating the magnitude of the phonon scattering rate and thermal conductivity.28,49,50Fig. 3c shows that the absolute value of the γ parameter of AC′ is close to that of AA′ in the range of 2–10 THz, indicating a weak influence of γ on the thermal conductivity of AC′ in that range. Similar to WP3, the total scattering rate of AC′ is higher than that of AA′ in the range from 2 to 10 THz as shown in Fig. 3d, which is consistent with the lower κ of AC′. The frequency differences of transverse optical (TO) and longitudinal optical (LO) modes also lead to the higher phonon scattering rate of AC′ as shown in Fig. S4 (ESI†). The atomic vibration amplitude of AC′ is also larger than that of AA′ indicating a higher lattice anharmonicity as shown in Fig. S5 (ESI†).51 The phonon mean free path (MFP) is an important reference for the modulation of thermal transport by nanostructuring.52 We calculate the phonon MFP based on the results of phonon relaxation time and group velocity.53 The MFP of AA′ is close to that of AC′ with the value around 300 nm as shown in Fig. S6 (ESI†). To modulate the thermal transport via nanostructuring, the scale of bilayer h-BN sample should be smaller than 1 μm. Overall, the variation of weak van der Waals interaction induces around 7% decrease of the thermal conductivity from AA′ (428 W m−1 K−1) to AC′ (397 W m−1 K−1). The lower κ indicates the AC′ configuration still exhibits acceptable heat dissipation efficiency if applied in the actual devices.
Based on above discussion, the softer ZO′ modes of AC increase the three-phonon scattering possibility and enhance the phonon scattering rate, resulting in a lower thermal conductivity. The soft phonon vibrations are generally led by a delocalized charge distribution or overlap of electronic wavefunction.54 Therefore, we next explore the possible reason for the soft ZO′ modes in AC′ by studying the electronic structures of AA′ and AC′. From the band structure of AA′ and AC′ in Fig. 4a and c, we can see that the conduction band of AC′ at the K point is obviously splitting unlike the 2-fold degenerate one of AA′. The projected density of states (PDOS) for the pz orbital of B atoms is significantly higher than that of N atoms in the conduction band as shown in the right panel of Fig. 4a and c. As can be seen, the conduction (valence) band of both AA′ and AC′ is contributed by B (N) atoms due to the stronger electronegativity of N compared with B. Hence, the stronger interaction between B atoms in different layers leads to a splitting in the conduction band of AC′. Furthermore, the pz states of B atoms in AC′ in the conduction band shift down compared with AA′ indicating a larger overlap of wavefunction for pz orbitals. The band gap of AC′ (4.21 eV) is also lower than that of AA (4.49 eV). The charge density between the two layers in AC′ is more delocalized than that of AA′ as shown in insets of Fig. 4, which is consistent with the larger overlap of wavefunctions.
Due to the decreasing of the bandgap, the Seebeck coefficient (S) and electrical conductivity (σ) of AC′ are also different from that of AA′ as shown in Fig. 4b and d. The absolute Seebeck coefficient value (|S|) of AC′ is lower than that of AA′ and the electrical conductivity is higher especially for the p-type, which is consistent with the smaller bandgap of AC′ and larger electron relaxation time as shown in Fig. S7a (ESI†). As a result, the power factor of AA′ is smaller than that of AC′ under p-type doping but larger under n-type doping as shown in Fig. S7b (ESI†). All the band splitting, downshift of PDOS and the delocalized charge distribution reveal the reason for the soft ZO′ mode in AC′ and further explain the lower thermal conductivity of the AC′-staggered bilayer h-BN than AA′ and AB (see Fig. S8 for AB, ESI†).
Fig. 5a shows the phonon band structure of AC′ with/without strains. As can be seen, the ZO′ frequency keeps increasing with the increase of the strain. That is to say, the stronger interlayer interaction by strain induces the higher frequency of ZO′ modes. Interestingly, unlike the monotonic change of the phonon frequency, the thermal conductivity of AC′ firstly decreases, then increases to the value even higher than that without strain as shown in Fig. 5b. The κ results based on the relaxation time approximation (RTA) have a similar trend to the iteratively solved data, but lower absolute values. It is because the RTA approach treats both normal and Umklapp scattering as the resistance for thermal transport, but the momentum conserving normal scattering does not resist the thermal transport.56,57 Hence, the RTA approach underestimates the thermal conductivity especially in materials with strong normal scattering and high thermal conductivity, such as h-BN and graphene. A large difference between the RTA and iterative way is significant at 0% and 18%, which indicate that normal scattering is more dominating in phonon scattering processes without strain or under a very large compressive strain.
To clarify the dominant reason for the thermal conductivity under various strains, we check the intermediate physical quantities for AC′ with strains. First, the WP3 value decreases monotonically with the strain increasing as shown in Fig. 5c, which is reasonably consistent with the increase of the phonon frequency. From eqn (S2) and (S3) (ESI†),48 the WP3 value just depends on the phonon dispersion relation. Therefore, the decreasing WP3 does not contribute to the decreasing thermal conductivity when the strain increases from 0% to 6%. Under a higher strain larger than 6%, the WP3 contribute to the increase of thermal conductivity to a higher extent. The anharmonic scattering rate (ASR) exhibits a similar variation trend to the thermal conductivity unlike WP3 (Fig. 5d). The ASR includes the effect in interatomic force constant as shown in eqn (S4) and (S5) (ESI†).33 The consistency between thermal conductivity and the ASR indicates that the anharmonic IFC or the phonon anharmonicity dominates the thermal conductivity more under the applied strain.
The −6% strain significantly enhances the anharmonic scattering of out-of-plane phonon modes (ZA and ZO′), while the further increasing strains larger than −6% gradually suppress the anharmonic scattering. The local maximum κ under the strain of 10% is also induced by the local minimum ASR as shown in Fig. S9 (ESI†). The final total scattering rate (Fig. 5e) has the same trend as ASR. The magnitude of γ parameter decreases with the increasing strain as shown in Fig. S10 (ESI†). The small γ under −18% further verifies the weak phonon anharmonicity consistent with the lower ASR value. For AA′ and AB configurations, there are obvious imaginary phonon frequencies near Γ point under a high compressive strain as shown in Fig. S11 (ESI†). We only investigate the κ of AA′ and AB with compressive strain below 6%. Both the κ of AA′ and AB exhibit a decreasing trend when the strain increases from 0 to 6% in Fig. S12 (ESI†). The ASR of AA′ reach the highest value under the 6% strain similar to the case of AC′ consistent with the smaller κ value. All the N (B) atoms of AA′ are in the top site inducing a stronger B–N van der Waals interaction, a higher ASR and a lower κ value under strain. Overall, the researchers should avoid any out-of-plane compressive strain on bilayer h-BN to keep a high heat dissipation efficiency for the devices.
To understand the layer effect on the thermal conductivity of sliding h-BN, we compare the intermediate quantities for AC′- and AA′-type tri-layer and bulk structure. Fig. 6 shows obvious differences in WP3 and the total scattering rate between AC′- and AA′-type tri-layers and bulk h-BN. The WP3 of AC′ is clearly larger than that of AA′ over the range of all phonon frequencies (Fig. 6a and d). The total scattering rate of AC′ is larger than that of AA′ mainly within 0–10 THz (Fig. 6b and e). More layers will effectively enhance the phonon scattering rate for AC′ with B–B patterns compared with AC′ with B–N patterns, thereby reducing the thermal conductivity significantly compared with the bilayer case. The MFP of the bulk structure is also obviously smaller than the that of the tri-layer structure from the cumulative κ in (Fig. 6c and f). For the tri-layer, the κ of AC′ is 13% lower than that of AA′ while the bilayer has only a 7% lower κ. More layers also enhance the phonon anharmonicity as shown by the change of ASR and γ in Fig. S15 (ESI†). In a short summary, κ will decrease with the increase in the number of layers of h-BN. The AC′-type configuration (B–B) always has the lowest κ in multilayer and bulk h-BN. From bilayer to bulk, the difference of κ between AA′ and AB also gets larger but is less obvious compared with AA′ and AC′ as discussed in the ESI† (Fig. S16, ESI†).
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3tc01531a |
This journal is © The Royal Society of Chemistry 2023 |