In-plane thermal conductivity of hexagonal boron nitride from 2D to 3D

Hexagonal boron nitride (h-BN), a two-dimensional layered material akin to graphene, has recently garnered significant attention due to its exceptional electronic, thermal, optical, and mechanical properties.^1–19 In particular, the weak van der Waals interactions allow a large freedom in precisely manipulating its thickness and thereby properties. For instance, the thickness dependence of bandgap in h-BN, transitioning from a 6.47 eV direct bandgap in monolayer form to a 5.95 eV indirect bandgap in bulk,^12,13 opens up possibilities for its use as a tunnel barrier or dielectric in electronic devices,²⁰ as well as in ultraviolet light emitters.^21,22 The optical properties of h-BN, such as second-order nonlinearities and hyperbolic dispersion, also largely depend on its thickness, offering great potential in the applications of nonlinear optics and photonics, especially in the mid-infrared range.^15,16 In the majority of relevant applications, thermal conductivity emerges as a pivotal factor influencing their heat dissipation and performance efficacy. Moreover, h-BN thin films have been actively considered as a promising heat spreader because of their high in-plane thermal conductivity and good dielectric characteristics.^23,24

Considering the varying thickness of h-BN in relevant applications, understanding the thickness dependence of its thermal conductivity is highly desirable. Up to now, very few studies have been devoted to investigating the thermal conductivity of h-BN with varying thickness, limited up to five layers.^1,5,6 The previous calculations based on the solution of the Peierls–Boltzmann transport equation (PBTE) reveal that, as the number of layers increases, its in-plane thermal conductivity decreases and quickly approaches the bulk value (∼443 Wm⁻¹K⁻¹) at around four layers.¹ A similar layer-dependence of thermal conductivity has been observed for graphene.² As the thickness further increases beyond five layers, how its in-plane thermal conductivity varies and recovers the bulk value remains elusive. Several experiments were performed to measure the in-plane thermal conductivity of h-BN samples with a thickness up to 20 nm,^7,25 which are all notably lower than the bulk values. In particular, the measurements by Jo et al.⁶ show that the 5-layer h-BN has a lower thermal conductivity than the 11-layer one. These results indicate that the thermal conductivity of h-BN may decrease below the bulk value at a nanoscale thickness, which, however, requires further exploration considering many factors in experiments (such as defects and contamination for the samples) can largely reduce its thermal conductivity. The transition of in-plane thermal conductivity from a monolayer to a bulk form is nontrivial. This is evidenced by a non-monotonic thickness dependence of thermal conductivity observed in other layered materials, such as MoS₂,²⁶ WTe₂,²⁷ black phosphorous,²⁸ Bi₂Te₃,²⁹ TaSe₂,³⁰ and TiS₃.³¹ This phenomenon is largely attributed to the interplay between phonon–phonon scattering and phonon-boundary scattering.

To comprehend thermal transport in layered materials like h-BN, atomistic simulations are often employed. The two common atomic-scale modeling methods include ab initio calculations grounded in density functional theory (DFT)³² and molecular dynamics (MD)^33,34 simulations based on empirical potentials.³⁵ DFT calculations, while offering high precision in describing interatomic interactions, are computationally intensive, thus constraining simulations to a few hundred atoms over picosecond durations. This computational expense makes it challenging, if not unfeasible, to simulate h-BN with far more than five layers using DFT calculations. On the other hand, MD simulations combined with empirical potentials, although generally computationally affordable, tend to be less accurate and have lower transferability compared to DFT calculations.

Recently, there has been a surge in developing accurate and efficient interatomic potentials using machine learning. Machine learning interatomic potentials (MLIPs), such as the neuroevolution potential (NEP),³⁶ moment tensor potential (MTP),³⁷ Gaussian approximation potential (GAP),³⁸ and deep potential (DP),³⁹ have demonstrated their effectiveness in various applications. Specifically, these MLIPs have been used in studying thermal conductivity and phonon transport in materials such as diamond,⁴⁰ crystalline silicon with vacancies,⁴¹ cubic boron arsenide,^42,43 wurtzite boron arsenide,⁴⁴ CoSb₃,^45, β-Ga₂O₃,⁴⁶ C₆₀,⁴⁷ graphene/borophene heterostructures,⁴⁸ bulk h-BN,⁴⁹ V^A–N binary compounds,⁵⁰ typical two-dimensional materials,⁵¹ graphene grain boundaries,⁵² and twisted van der Waals structures,⁵³ offering accuracy comparable to DFT but with significantly lower computational costs.

In this study, we investigated the in-plane thermal conductivity of h-BN from its two-dimensional (2D) to three-dimensional (3D) forms using the MLIP developed within the NEP framework. Our NEP model demonstrated a high precision comparable to DFT in predicting the lattice thermal conductivity and phonon properties across various h-BN configurations, including monolayer, multilayer, and bulk. By utilizing the homogeneous nonequilibrium molecular dynamics (HNEMD)^54,55 method combined with the NEP model, we calculated the in-plane thermal conductivity of h-BN with a thickness up to approximately 100 nm (300 layers). We observed a rapid reduction in its thermal conductivity from a monolayer to a multilayer configuration, eventually stabilizing at its bulk value after exceeding four layers. This phenomenon was further explained using the thickness dependence of phonon group velocity and lifetime determined by the spectral energy density (SED)⁵⁶ method. In particular, the weak thickness dependence of phonon lifetime in nanoscale h-BN thin films can be attributed to its exceptionally high phonon focusing along the in-plane direction. This research fills the gap in understanding the evolution of thermal conductivity from 2D to 3D h-BN, which is helpful for the thermal design and performance optimization of h-BN-based electronics and optoelectronics.

The MLIP within the NEP framework was originally proposed by Fan et al. as a promising tool to study thermal transport with high accuracy and low cost.³⁶ It has been improved later and the version we used in this paper is the NEP4 model.^57,58

The NEP approach employs a single-layered feedforward neural network for predicting the site energy of atoms in a material. A key feature of NEP is its implementation of unique descriptors for every atom, which combine radial and angular components for comprehensive representations of atomic environments. A NEP model is trained via the separable natural evolution strategy,⁵⁹ which focuses on reducing a loss function. This loss function is formulated as a cumulative weighted sum that includes the root mean square errors (RMSEs) of energy, force, and virial coefficients, combined with components that provide L₁ and L₂ regularization. The “virial” here measures the mechanical stress at an atomic scale, defined as $W i= ∑ j ≠ i r i j ⊗ ∂ U j ∂ r j i$⁠, where i and j denotes the atom index, $U j$ represents the potential energy of atom j, and $r i j≡ r j− r i$ is the vector denoting the distance between atoms j and i. More details about the NEP methodology can be found in previous work,³⁶ and the definitions of all hyperparameters are provided in Table SI in the supplementary material.

In this study, the training database was constructed by configurations of monolayer, 2-layer, 3-layer, 4-layer, and bulk h-BN, informed by evidence that its thermal conductivity decreases as layer count increases and stabilizes at the bulk level by around four layers.¹ Specifically, the training database is composed of 230 bulk h-BN, and 20 each for monolayer, 2-layer, 3-layer, and 4-layer h-BN. These configurations were selected from ab initio molecular dynamics (AIMD) simulations at 300 K using a 1 × 1 × 1 k-point grid, capturing snapshots at intervals of every 50 time steps. Each configuration consists of 200 atoms for bulk h-BN, 50 atoms for monolayer h-BN, 100 atoms for 2-layer h-BN, 150 atoms for 3-layer h-BN, and 200 atoms for 4-layer h-BN, giving a total of 56 000 environments. These configurations were further refined through self-consistent field calculations with a 2 × 2 × 2 k-point grid to obtain the total energy, atomic forces, and virial. Leveraging this methodology, we constructed a testing database containing 20 configurations each for monolayer, 2-layer, 3-layer, 4-layer, and bulk h-BN. All DFT calculations utilized the projector augmented wave (PAW)⁶⁰ method within the Vienna Ab initio Simulation Package (VASP),⁶¹ employing the local density approximation (LDA)⁶² for the exchange-correlation functional and a cutoff energy of 600 eV.

After determining the training and testing databases, we employed the GPUMD package⁵⁷ to train the NEP model for h-BN. The primary hyperparameters for NEP training are outlined in Table I. The NEP model was trained over 100 000 generations, with loss functions for each component beginning to show signs of convergence around 50 000 generations and achieving full convergence by the end of the training process [Fig. S1 in the supplementary material].

For all simulations employing the HNEMD approach, the system is first equilibrated for 0.5 ns at 300 K based on the NVT (keeping the number of atoms, volume, and temperature constant) ensemble controlled by the Nose–Hoover chain thermostat,⁶³ with a time step of 1 fs. Subsequently, the HNEMD method is applied to calculate the thermal conductivity of h-BN within the NVT ensemble, employing a time step of 1 fs and a temperature maintained at 300 K. The driving force parameter was set at F_e = 0.1 μm⁻¹ for h-BN across monolayer, multilayer, and bulk forms. As the number of h-BN layers increases, the size of the simulation cells we use also grows, whereas the production time decreases. Specifically, for h-BN with one to four layers and the bulk, the production time was chosen at 3 ns. For multilayer h-BN structures with 10–100 layers, the production time was reduced to 1 ns, and further reduced to 0.6 ns for 300-layer h-BN. For all configurations of h-BN, we use simulation cells containing over 20 000 atoms to ensure that our simulations are conducted on a sufficiently large scale. Six separate HNEMD simulations were conducted to obtain an average thermal conductivity, with an estimation of statistical error.

We calculate the thermal conductivity by resolving the linearized PBTE, as implemented in SHENGBTE.⁶⁵ Detailed information of PBTE calculations is available in previous works.^43,49,66

To obtain the thermal conductivity of h-BN, both second-order and third-order interatomic force constants (IFCs) were first calculated using the NEP and DFT via the finite displacement approach. A 5 × 5 × 1 (5 × 5 × 2) supercell was used for the monolayer and multilayer h-BN (bulk h-BN). The IFCs determined by the NEP were executed using the LAMMPS package.⁶⁷ The DFT calculations of IFCs utilized a 3 × 3 × 1 (3 × 3 × 3) k-mesh for monolayer and multilayer h-BN (bulk h-BN), with the cutoff energy being 600 eV. Convergence criteria for the energy and force in DFT calculations were set at 10⁻⁸ eV and 10⁻⁶ eV/Å, respectively. We employed the HiPhive package⁶⁸ to enforce the rotational invariance conditions on the second-order IFCs of monolayer and multilayer h-BN, which is crucial for accurately capturing the quadratic phonon dispersion of ZA branch near the Gamma point. To ensure thermal conductivity convergence, the third-order IFC calculations accounted for up to the 5th (12th) nearest neighbors for monolayer h-BN (multilayer and bulk h-BN). Using the derived IFCs, we calculated the lattice thermal conductivities of h-BN by iteratively solving the linearized PBTE using SHENGBTE.⁶⁵ This process utilized different q-mesh sizes to ensure convergence, with 120 × 120 × 1 for monolayer, 50 × 50 × 1 for multilayer, and 24 × 24 × 9 for bulk.

We first evaluated the accuracy of the developed NEP model using the RMSEs of atomic energy, force, and virial for h-BN with varying layer numbers. Figure 1 reveals that the RMSE for monolayer, 2-layer, 4-layer, and bulk h-BN are comparable. For the training database, the highest RMSE values are 0.21 meV/atom for atomic energy, 24.18 meV/Å for atomic force, and 4.12 meV/atom for virial. In the testing database, these RMSE values slightly increase to 0.24 meV/atom for atomic energy, 26.49 meV/Å for atomic force, and 4.53 meV/atom for virial. We note that these RMSE values for atomic forces are relatively lower than those reported NEP models for other layered materials, such as graphene (46.72 meV/Å),⁷⁰ MoS₂ (53 meV/Å),⁵³ and PbTe (38 meV/Å).⁷¹ 

We continue to evaluate the accuracy of the NEP in predicting the phonon transport properties and thermal conductivity of h-BN. Figure 2 shows that the phonon dispersions predicted by NEP almost overlap with DFT results for monolayer, 2-layer, 4-layer, and bulk h-BN, indicating that the NEP can well predict the second-order IFCs. Furthermore, the NEP exhibits a DFT-level accuracy in predicting the in-plane thermal conductivities of all these structures [Fig. 3] and the cross-plane ones of bulk h-BN [Fig. S2 in the supplementary material] within the temperature range from 100 to 1000 K, demonstrating its capability to reproduce third-order IFCs. All these results demonstrate the high precision of the NEP model in reproducing the ab initio potential energy surfaces of h-BN across monolayer, multilayer, and bulk forms, a capability critical to exploring the thickness-dependent thermal conductivity of h-BN.

Employing the developed NEP, we next look into how thickness influences the in-plane thermal conductivity of h-BN. Figure 4 shows the results calculated by HNEMD combined with the NEP model for a thickness up to 100 nm at 300 K, in comparison with the reported theoretical and experimental results, as well as our PBTE results. As shown in Fig. 4, as the thickness increases, the predicted thermal conductivity rapidly decreases and approximates the bulk value (496.7 ± 29.7 W m⁻¹K⁻¹) at four layers. This observation is in agreement with both our PBTE results and the reported PBTE¹ and MD results.⁵ As the thickness further increases up to 100 nm, the thermal conductivity remains almost unchanged. This behavior contradicts the general expectation that phonon-boundary scattering can strongly suppress the thermal conductivity of nanoscale thin films, which has been commonly observed in thin films such as MoS₂,²⁶ WTe₂,²⁷ and TiS₃.³¹ 

We notice that our thermal conductivity calculated by PBTE exhibits a similar thickness dependence but deviates significantly from the values predicted by HNEMD. The PBTE value of monolayer h-BN calculated considering only third-order anharmonicity is significantly higher than the HNEMD result, while for multilayer and bulk h-BN, it is lower than the HNEMD result. This discrepancy could be attributed to the effect of both four-phonon scattering and temperature dependence of IFCs on the thermal conductivity. According to previous works, four-phonon scattering reduces thermal conductivity,^76,77 whereas considering the temperature dependence of IFCs tends to increase it.⁷⁸ When using temperature-dependent IFCs, four-phonon scattering still significantly reduces the thermal conductivity of monolayer graphene by around 50%.⁷⁹ For monolayer h-BN, which has a similar structure as monolayer graphene, further considering both four-phonon scattering and the temperature dependence of IFCs may also result in an overall reduction in thermal conductivity. In contrast, the impact of four-phonon scattering on thermal conductivity is much weaker in multilayer graphene and graphite due to the breaking of reflection symmetry selection rules and the splitting of the flexural acoustic branch.⁷⁶ As reported by Li et al.,⁸⁰ four-phonon scattering results in a small reduction of ∼15% in the cross-plane thermal conductivity of graphite. Similarly, our calculation⁴⁹ shows that further including four-phonon scattering decreases the thermal conductivity of bulk h-BN at 300 K only by 6.91% and 3.50% along in-plane and cross-plane directions, respectively. Therefore, the underestimation of thermal conductivity in multilayer and bulk h-BN predicted by PBTE compared to the HNEMD results could be attributed to the temperature dependence of IFCs. As reported by Yang et al., when both three-phonon and four-phonon scattering are included, considering the temperature dependence of IFCs results in ∼50% enhancement in the thermal conductivity of UO₂ at 300 K.⁷⁸ The research by Gupta et al.⁸¹ also indicates that considering the temperature dependence of IFCs is vital to understanding and predicting the thermal conductivity of bulk MoSe₂ and WSe₂. The deviations between the PBTE and HNEMD results highlight the significance of machine learning interatomic potentials, allowing MD simulations with a DFT-level accuracy. MD simulations naturally include all-order anharmonicity and temperature effects, whose calculations are expensive or even prohibitive in the PBTE approach. In particular, MD simulations combined with machine learning interatomic potentials have demonstrated their success in understanding the effect of strong anharmonicity on lattice dynamics and thermal transport.^82,83

To understand the thickness dependence of thermal conductivity, we calculated the spectral thermal conductivity of h-BN across monolayer, multilayer, and bulk at 300 K. As shown in Fig. 5, the spectral thermal conductivity decreases significantly from monolayer to 2-layer, mainly below 10 THz, while it varies little as the thickness further increases. The phonon modes below 10 THz are mainly acoustic modes [see Fig. 2], among which the flexural phonons (ZA) contribute most to the overall thermal conductivity, as predicted by Lindsay and Broido.¹ In Fig. 6(a), we compare the thermal conductivity contribution from the acoustic branches for the monolayer and 2-layer h-BN, which are calculated by the NEP-based PBTE approach. Clearly, the suppression of thermal conductivity in the 2-layer h-BN as compared to the monolayer one primarily results from the ZA modes. The main reason for the lower thermal conductivity of the ZA branch in two-layer h-BN is due to its higher scattering rates [Fig. 6(b)] and the reduced number of ZA phonons [Fig. 6(c)]. In contrast, the group velocities of ZA phonons in 2-layer h-BN are very close to those in the monolayer one [Fig. 6(d)].

We continue to look into the phonon transport details of h-BN with a larger thickness. As shown in Fig. 7(a), the phonon group velocities of acoustic phonons change little from 2-layer to bulk, indicating its negligible effect on thermal conductivity. Figure 7(b) further shows the phonon lifetime calculated by the SED method for monolayer, 2-layer, 10-layer, 50-layer, and bulk h-BN at 300 K. The overall shape of phonon dispersion calculated by SED is apparent for each structure [Fig. S3 in the supplementary material], which highlights the precision of phonon lifetime determined by fitting the SED curve. One can see that the phonon lifetime of monolayer h-BN is slightly larger than others through the entire frequency range while the increase of thickness above two layers results in negligible suppression in the phonon lifetime.

The little suppressed phonon lifetime of nanoscale thin films compared to its bulk indicates weak phonon-boundary scattering in h-BN, which can be attributed to the extreme anisotropy of phonon transport. Figure 8 illustrates the ratio of group velocity in the z-direction (cross-plane) to that in the y-direction (in-plane) v_z/v_y for bulk h-BN. We find that v_z is significantly smaller compared to v_y, with only 12.17% of phonon modes having a v_z/v_y ratio greater than 1, indicating strong phonon focusing along the in-plane direction. Based on our PBTE results, the phonon modes with v_z/v_y > 1 account for only 0.29% of the in-plane thermal conductivity in bulk h-BN. Even for phonon modes with v_z/v_y > 0.1, the corresponding contribution is merely 11.97%.

A further demonstration was performed to calculate the in-plane thermal conductivity of h-BN thin films by including the phonon-boundary scattering in the PBTE calculation. To enable a direct comparison with the thermal conductivity value acquired through HNEMD (496.7 ± 29.7 W m⁻¹K⁻¹), we adjusted the value obtained by PBTE from 276.5 to 499.4 W m⁻¹K⁻¹ by scaling the phonon scattering rate. We consider the thin films with a thickness down to 10 nm because the phonon dispersion of an even thinner film can differ significantly from the bulk case. Figure 4 presents the thermal conductivity obtained with varying values of specularity parameter p for the phonon-boundary scattering. Because the free-boundary has been used in the MD simulations, the value of p is expected to be close to 1, which corresponds to total specular reflection at the boundary and causes no reduction in thermal conductivity. For p = 0.9, the thermal conductivity exhibits a very weak thickness dependence, with a small reduction of 4.39% as the thickness decreases to 10 nm. Adjusting p to 0.5 results in a slightly larger reduction (14.24%) in the thermal conductivity for a h-BN film with a thickness of 10 nm compared to the bulk one. Even for p = 0, corresponding to the case of complete diffusion scattering, the thermal conductivity of a 10-nm thick h-BN film is still as high as 371.1 W m⁻¹K⁻¹, which is only 25.70% lower than the bulk value (499.4 W m⁻¹K⁻¹).

In summary, we have developed a high-quality MLIP within the NEP framework for h-BN. This NEP model achieves a DFT-level precision in predicting lattice thermal conductivity and phonon properties across different thicknesses of h-BN from a monolayer to its bulk material. By combining MD simulations and the NEP model, we predicted the in-plane thermal conductivity of h-BN with a thickness from monolayer to approximately 100 nm (300 layers), finding that the thermal conductivity decreases rapidly from a monolayer and stabilizes at its bulk value beyond four layers. This convergence of in-plane thermal conductivity is due to the negligible changes in both phonon group velocity and lifetime as thickness grows above four layers. The weak thickness dependence of phonon lifetime in its nanoscale thin films originates from its significant in-plane phonon focusing. In particular, phonon modes with a velocity ratio of v_z/v_y > 1 account for only 0.29% of the in-plane conductivity in bulk h-BN.

Our research reveals the evolution of in-plane thermal conductivity in h-BN from its monolayer to bulk, which deepens the understanding of phonon transport from 2D to 3D and will benefit the thermal design and performance optimization of h-BN-based electronics and optoelectronics.

See the supplementary material for figures depicting the evolution of loss function components during training, the cross-plane thermal conductivity of bulk h-BN calculated by the NEP and DFT, the SED spectra calculated by MD simulations based on the NEP for h-BN at 300 K, and definitions of hyperparameters for the NEP training.

We acknowledge support from the Excellent Young Scientists Fund (Overseas) of Shandong Province (Grant No. 2022HWYQ-091), the Taishan Scholars Program of Shandong Province, the Natural Science Foundation of Shandong Province (Grant No. ZR2022MA011), and the Initiative Research Fund of Shandong Institute of Advanced Technology.

The authors have no conflicts to disclose.

Jialin Tang: Data curation (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead). Jiongzhi Zheng: Methodology (supporting); Writing – review & editing (supporting). Xiaohan Song: Methodology (supporting); Writing – review & editing (supporting). Lin Cheng: Supervision (equal); Writing – review & editing (equal). Ruiqiang Guo: Conceptualization (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.