The escalating demand for nanoscale thermal management in advanced electronic devices exhibits the importance of interfacial heat transfer, which critically governs device performance and reliability [1, 2]. Interfacial thermal conductance (ITC), a key parameter characterizing heat transport across interfaces, can be tuned by structural configurations, chemical composition, and defect [3, 4]. For instance, electron-phonon coupling dominates heat transfer at metal-insulator interfaces [5, 6], while wettability and interfacial morphology dictate solid-liquid thermal transport [7].

In recent years, two-dimensional lateral heterostructures have demonstrated unique properties and potential applications to tunnel field-effect transistors, memory devices, and optoelectronic devices [8–10]. Due to the polymorphism of borophene, different phases of the structure can coexist in experimental preparations. As a result, both theoretical and experimental investigations of borophene-based two-dimensional lateral heterostructures have garnered increasing attention [11, 12]. Studies have shown that distinct borophene phases exhibit different phonon dispersion relations, leading to variations in their thermal transport properties [11, 13]. The impact of planar interfaces in borophene on thermal transport is multifaceted: interfaces between different borophene phases can induce a phenomenon of thermal rectification, i.e., a directional asymmetry in thermal transport, thereby providing opportunities for the development of novel thermal management devices [14]. Additionally, the discontinuity of atomic structures at the interfaces generates thermal resistance, impeding heat transfer [15]. The extent of this thermal resistance is closely related to the specific interface structure. Taking graphene-borophene heterostructures as an example, phase boundaries significantly influence phonon transport at interface, thus modulating the system’s ITC [16].

Concurrently, the synthesis of borophene has remained a challenging issue in laboratory settings. The difficulty in synthesizing borophene arises partly from the ability of boron atoms to form two- or three-center B-B bonds, which are more prone to interaction and the formation of polyhedral structures, making the bonding patterns more complex than those between carbon atoms. This often results in various defects during the synthesis process [17–19]. The random vacancy defects in borophene sheet enhance phonon inelastic scattering, thereby reducing thermal conductivity [20]. Researches indicate that interface engineering is an effective means of regulating ITC [21]. For example, the introduction of interface coupling agents in composite materials can enhance thermal conductivity [22], while optimizing the interface structure represents another effective method to improve thermal conductivity [23]. Commonly, the incorporation of defects has been considered to significantly reduce in-plane thermal conductivity [24–26]. However, recent works have indicated that defects can also promote phonon transmission at the asphalt/graphene interface, leading to a substantial increase in ITC [27]. In MOF/TiO₂ heterostructures, hydrogen bond channels is regarded as an important factor for facilitating interfacial heat transfer [28]. Furthermore, certain defects can also guide phonons along specific paths to suppress disordered scattering at the graphene/h-BN interface [29]. Therefore, the influence of the interface on the thermal transport of the low-dimension materials is complexed.

Despite extensive research on various lateral heterostructures, the literature on borophene-based heterostructures (BHs) remains relatively sparse. The polymorphism of borophene in these heterostructures is often overlooked. Therefore, a comprehensive understanding of the role of borophene polymorphism in lateral heterostructures is crucial. Based on prior studies, this work employs molecular dynamics (MD) simulations utilizing a machine-learning-based moment tensor potential (MTP) to investigate the thermal transport properties of the β ₁₂/χ ₃ mixed-phase borophene interface and the influence of defects on the ITC. Through the analysis of interfacial thermal conductance, spectral heat flux, and stress distribution characteristics, this study aims to enhance the understanding of the mechanisms governing interfacial thermal transport and provide theoretical guidance for the design of high-performance thermal management materials.

The structural optimization was performed using density functional theory (DFT) calculations with the Vienna Ab-initio Simulation Package (VASP). The Generalized Gradient Approximation (GGA) was employed, utilizing the Perdew-Burke-Ernzerhof (PBE) functional to describe electron-electron interactions. The cutoff energy for the plane-wave basis set was set to 500 eV. The convergence criteria for energy and forces were set to 10^–6 eV/Å and 10^–3 eV/Å, respectively. A vacuum layer with a thickness of 20 Å was applied in the direction perpendicular to the borophene sheet. The lattice constants of β ₁₂ (χ ₃) borophene along the zigzag and armchair directions are 8.41 (5.06) Å and 2.90 (2.93) Å, respectively, which are in excellent agreement with previous studies [30,31]. The atomic configuration of the two-dimensional β ₁₂/χ ₃ borophene heterostructures are shown in Figure 1A. Given that the lattice constants of the two borophenes along the armchair direction are nearly identical, we selected the zigzag direction perpendicular to it as the interface normal. The residual stress due to lattice mismatch in this case is less than 0.35 GPa. Figures 1B,C display the single-vacancy (SV) and double-vacancy (DV) defects introduced into the original BHs, hereafter referred to as SV-BHs and DV-BHs, respectively. These defects were introduced by removing one boron atom and two adjacent boron atoms, respectively, from the interface of the system, with the missing boron atoms highlighted in purple. The length of the BHs considered in the x-direction is 20.2 nm, and in the y-direction, it is 5.8 nm.

All MD simulations in this work were performed using the LAMMPS molecular dynamics simulator [32], which has been widely utilized for classical MD simulations. A machine learning-based interatomic potential model, the moment tensor potentials (MTPs) [33], based on the generalized linear framework, was employed to describe the atomic interactions within the B/B system. We have constructed an Ab-initio molecular dynamics (AIMD) dataset consisting of approximately 3,000 configurations for the development of the MTP. Supplementary Figure S1 presents a comparison of the phonon dispersion calculated via DFT and MTP for the χ ₃ and β ₁₂ structures. The results demonstrate that the phonon dispersion obtained from MTP calculations is in excellent agreement with the DFT results, thereby validating the accuracy of our machine learning model. In the non-equilibrium MD (NEMD) simulations, a schematic of the simulation domain is shown in Figure 1A. To prevent the system from undergoing translational motion during the MD simulation, atoms within the gray regions at both ends, approximately 5 Å in length, were fixed. The time step for the entire simulation process was set to 0.5 fs, and periodic boundary conditions were applied in all three directions. To establish a linear temperature distribution within the system, atoms in the red and blue regions were respectively subjected to heat baths at temperatures of T₀ + ΔT and T₀ - ΔT, where T₀ and ΔT were set to 300 K and 25 K, respectively.

Prior to conducting the thermal transport simulations, the initial borophene heterostructure was optimized through energy minimization using the conjugate gradient method. The atomic velocity distribution was initialized to a Gaussian distribution corresponding to 50 K. Subsequently, the temperature was gradually increased from 50 K to 300 K under zero pressure via an NPT relaxation process. The system was then relaxed under NPT ensemble conditions at 300 K and zero pressure to minimize structural energy. Finally, an NVE ensemble was applied to all atoms (except those that were fixed), and the steady-state heat flux was obtained over a 5 ns period. The temperature of atoms at any given moment within the system is calculated using the following Equation 1: T = m ∑ i = 1 N V i 2 3 N i k B (1)

In the equation, T represents the atomic temperature; m denotes the atomic mass; V _i is the velocity of the i-th atom; N _i is the total number of atoms; and k _B is the Boltzmann constant. Figure 2A presents the temperature distribution and energy over the final 5 ns, where the linear energy curve confirms the steady-state heat transfer. The slopes of these curves were used to compute the steady-state heat flux (J). The model was divided into 50 segments along the length direction, and the average temperature of atoms within each segment was calculated. The resulting temperature distribution along the length direction of the model at 300 K is shown in Figure 2B. Fourier’s law of heat conduction was employed to calculate the interfacial thermal conductance of the monolayer borophene heterostructures with different configurations by Equation 2: κ = J A · ∆ T (2)

Here, J denotes the heat flux along the free cross-section, and A represents the cross-sectional area through which the heat flux passes. In the present study, the interlayer distance of the monolayer borophene sheets is fixed at 4.26 Å ΔT indicates the temperature gradient across the heterostructure.

The finite size of the simulation cell inherently constrains the maximum attainable phonon mean free path (MFP), necessitating the consideration of finite size effects. Figure 3A illustrates the dependence of the interfacial thermal conductance of the pristine BHs on their length at 300 K. With the width fixed at 5.8 nm, the ITC of BHs with lengths of 20, 30, 40, 50, and 60 nm was calculated. As the length increases, the ITC of the BHs rises from 6.56 GWK⁻¹m^-2 to 18.31 GWK⁻¹m^-2, exhibiting a relatively high interfacial thermal transport [34,35], which is consistent with the trend observed in the literature for the thermal conductivity of borophene nanoribbons [15]. When the sample size is small, especially when it is smaller than the phonon MFP, the ITC increases sharply with length. The relationship between them is described as follows [36]: 1 / κ L = 1 + ∧ / L / κ ∞ , where κ L denotes the lattice thermal conductivity at a finite length L, while κ ∞ represents the lattice thermal conductivity at infinite length. ∧ refers to the effective phonon MFP. This behavior is attributed to the enhanced contribution of low-frequency, long-wavelength phonons to thermal transport. These phonons are more easily transmitted through the interface and exhibit lower inelastic scattering, resulting in a higher ITC [37].

To quantitatively assess the impact of defects on the interface, the ITC of pristine and defective BHs were calculated and presented in Figure 3B. The results indicate that both SV and DV defects in borophene significantly reduce the ITC values. As defects involve atomic vacancies, phonon transport at defect sites leads to local phonon scattering, which introduces local thermal resistance. Moreover, BHs with SV defects exhibit a relatively higher ITC value of 3.14 GW K⁻¹m⁻² approximately double that of BHs with DV defects. The observed variations in ITC are attributed to the geometric discrepancies at the interface and the atomic interactions between the different materials of the heterostructure [38]. The following discussion will delve into the mechanisms by which defects influence the interface behavior.

It is important to note that the atomic interactions at the interface predominantly govern the phonon transport across the interface [39]. To further investigate the vibrational differences at the three interfaces, atomic vibration analysis was performed on atoms located within 2 nm of the interface in the frequency domain. The overlap of the vibrational density of states (VDOS) plays a crucial role in the thermal transport between the two contacting materials. The VDOS of boron atoms on either side of the interface in the borophene heterostructure was calculated using the Fourier transform of the velocity autocorrelation function (VACF) by Equation 3 [40,41]: DOS ω = 1 2 π ∫ 0 τ v j 0 v j t e i ω t d t (3)

Here, ω represents the frequency, v denotes the velocity of carbon atoms, and ⟨…⟩ indicates the time and atomic number averaged velocity autocorrelation function, with τ being the computation time duration. The ensemble average in Equation 3 is achieved by performing a time average of the system’s atomic velocities over a 250 ps time span, extracting the velocity every 10 fs from the simulation. Figure 4A presents the VDOS of β ₁₂ and χ ₃ borophene at the pristine BHs interface, with the selected frequency range (0–40 THz) encompassing the phonon vibrational range of the borophene-based material, primarily originating from the vibrations of covalent B-B bonds. The similar density of states in both regions and the broad frequency range facilitate significant phonon excitation at the interface during thermal transport, thereby enhancing the thermal conductivity. Figures 4B,C illustrate the VDOS curves of the BHs with SV and DV defects. At certain peaks within the VDOS, splitting into multiple peaks occurs, indicating strong phonon scattering at the defects. To quantify the VDOS overlap, we calculated the overlap parameter (S) by Equation 4 as follows [34,42]: S = ∫ 0 ω max ω A ω d ω (4)

In this context, A(ω) represents the coherent part at frequency ω, which is displayed in the form of the S value in Figure 4. A higher overlap factor indicates a greater degree of phonon matching at the β ₁₂/χ ₃ interface. The results demonstrate that, compared to pristine borophene, the VDOS curve of borophene at the interface of SV-BHs shows a greater overlap in the low-frequency modes, within the vibrational frequency range of β ₁₂ and χ ₃ borophene being comparable, yielding an S value of 83.7%. In contrast, the variation for DV-BHs is smaller, with an S value of 74.8%, suggesting that the phonon transport channels at the interface of SV-BHs are more abundant than those in DV-BHs.

Since the phonon transport in borophene exhibits anisotropy across different directions, in-plane transport is significantly faster than out-of-plane transport [43]. The enhanced thermal conductivity between the β ₁₂/χ ₃ borophene heterostructures can be attributed to the higher overlap between the in-plane and out-of-plane VDOS curves of borophene. To better understand the effects of defects, the total VDOS curve of BHs, as shown in Figure 5, is separated into in-plane and out-of-plane modes. The atomic in-plane vibrational frequency range and intensity of BHs are notably larger than those of the out-of-plane vibrational modes, with the 18–40 THz frequency range predominantly corresponding to in-plane vibrations. Moreover, more characteristic peaks of in-plane vibrations are excited, which are primarily dominated by out-of-plane vibrations. Upon further discussion, the VDOS behavior of defective borophene differs. For defective borophene, as seen in Figures 5B,C, the overlap region of SV-BHs and DV-BHs VDOS does not exhibit significant changes. Simultaneously, the phonon matching degree at the interface under the SV defect is also higher, with an S value of 31.4%, exceeding the 29.5% observed under the DV defect. This agrees well with the general trend where better phonon matching typically results in higher ITC. Unfortunately, the S value for the pristine heterostructures is smaller than that for the defective structures, which can not explain the fact that it has higher ITC.

Further investigation was conducted using the spectral thermal flux (SHC) decomposition method to calculate the phonon transmission function, aiming to explore phonon heat transport. By analyzing the force-velocity correlation between atoms near the interface, the transmission function from material A to material B can be calculated by Equation 5 as follows [44,45]: Γ ω = 2 k B ∆ T M ∆ t s Re ∑ i ∈ A F → ̂ i ω · v → ̂ i ω * (5)

Here, k _B denotes the Boltzmann constant, M represents the number of samples, Δt _s is the sampling interval, and F _ij refers to the Fourier-transformed interatomic forces between atom-i in the material A region and atom-j in the material B region. The interatomic force and velocity data are recorded within a 2 nm region centered around the interface. In the calculations, we set Δt _s = 10 fs and M = 500000, with the results shown in Figure 6. The transmission curves exhibit similar shapes, reflecting the β ₁₂ and χ ₃ borophene in the heterogeneous structure. The transmission curves corresponding to BHs exhibit the highest peak, followed by SV-BHs, and then DV-BHs, consistent with the order observed in the ITC. Since phonons are the primary heat carriers, a higher transmission coefficient indicates that more phonons are able to traverse the interface at that frequency, meaning a greater amount of heat can be transferred through the interface. Meanwhile, a redshift is observed in DV-BHs. The defects lead to an increase in the atomic vibration frequency in the interface region (phonon hardening), which subsequently reduces the number of low-frequency phonons. The decreased contribution of low-frequency phonons diminishes thermal conductance [25]. According to Hooke’s law, the relationship between bond strength and phonon frequency can be expressed as ω = 1 2 π k μ , where ω is the phonon frequency, k is the force constant of the chemical bond, and μ is the reduced mass of the bonded atoms. Thus, the red shift of the peak position indicates a weakening of the bond strength at the interface, and weaker atomic bonds are detrimental to thermal conduction. This observation is consistent with the results shown in Supplementary Figure S2. We computed the stress-strain relationship near the heterogeneous structure interface under uniaxial tensile strain along the thermal transport direction. The strain is defined as ε = ΔL/L ₀, where ΔL represents the change of the length and L ₀ is the initial length. The stress monotonically increases with strain until it rapidly decreases upon reaching the tolerance limit. In the initial linear increasing segment of the stress-strain curve, the slope represents the Young’s modulus (E), which can be utilized to characterize the bond strength. In this context, within the heterogeneous structure, ITC exhibits a positive correlation with the bonding strength near the interface.

The improvement in density of states overlap facilitates thermal transfer, whereas excessive lattice distortion (quantified by von Mises stress) induces strong phonon localization effects. Therefore, to address the apparent contradiction between phonon compatibility and the degradation of thermal transport properties, we also calculated the von Mises stress distribution near the interface of three BHs (within 6 nm range). The atomic von Mises stress is defined by Equation 6 as follows [35]: σ = 1 2 σ x x − σ y y 2 + σ y y − σ z z 2 + σ x x − σ z z 2 + 6 σ x y 2 + σ y z 2 + σ x z 2 (6)

In this context, σ x x , σ y y and σ z z represent the normal stresses along the x, y, and z directions, respectively; σ x y , σ y z and σ x z denote the shear stresses. The atomic von Mises stress is employed to characterize the plastic behavior of an individual atom or a localized region. As shown in Figure 7, the stress distribution near the interface of the original BHs is relatively uniform, with excellent matching between the two materials in the heterogeneous structure. However, the existence of SV and DV defects result in a significant stress concentration near the interface. The degree of stress concentration in the DV-BHs near the interface is significantly higher than that in the SV-BHs. In the case of DV-BHs, after structural relaxation, the heterostructure adopts a chain-like connectivity, with a pronounced stress concentration. The von Mises stress of certain atoms at the interface reaches up to 40 GPa. Typically, higher stress concentrations inhibit atomic vibrations and increase phonon scattering at the interface, further deteriorating the thermal properties of the borophene-based heterogeneous structure. Consequently, the defects result in the reduction of ITC of β ₁₂/χ ₃ borophene heterostructures, and the DV-BHs has the lowest ITC.

In this study, we investigated the influence of defect on the interfacial thermal conductance of β ₁₂/χ ₃ borophene lateral heterostructures, employing first-principles calculations and Molecular Dynamics simulations based on MTP potential. The results demonstrate that the ITC of the pristine borophene heterostructures exhibits a pronounced size effect, significantly increasing with length (from 6.56 GW K⁻¹m^-2 to 18.31 GW K⁻¹m^-2). The introduction of single vacancies and double vacancies causes a reduction in ITC. The ITC of DV-BHs decreases to 1.57 GW K⁻¹m^-2. The phenomenon can not be explained by the overlap between the in-plane and out-of-plane VDOS at the interface. Subsequently, by analyzing spectral thermal flux and stress distribution, we revealed that the ITC is regulated by a competition between atomic bond strength and stress, which exacerbates phonon scattering. Our study provides theoretical guidance for the design of thermal management materials in nano-device applications.