Thermal boundary conductance (TBC) is important for heat dissipation in light-emitting diodes (LEDs). In this study, we predicted the TBC between the high thermal conductivity boron arsenide (BAs) and silicon (Si) by nonequilibrium molecular dynamics (MD) simulations. From the thermal conductivity accumulation function with respect to phonon frequency, the dominant phonon frequencies for heat conduction in BAs are extremely different from those in Si. However, our nonequilibrium MD simulations indicated that the TBC of the BAs/Si interface was still high compared to most other interfaces, even though there was a major frequency mismatch in the thermal conductivity accumulation function between BAs and Si. The primary reason for the high TBC is the overlap of phonon density of states between BAs and Si in the frequency range of 5–8 THz. The range of predicted TBC of the BAs/Si interface was between 200 and 300 MW/m^{2} K in the temperature range of 300–700 K, and the values of the TBC were not sensitive to the temperature. We also found that the TBCs in Si/BAs and Si/Ge interfaces were close to each other considering the simulation uncertainty. This work indicates BAs as an excellent material for heat dissipation across the interfaces.

## I. INTRODUCTION

The solid-state light-emitting diode (LED) is a promising technology to efficiently reduce electricity consumption. One of the major issues associated with LED is the heat generation and dissipation in the device.^{1} The accumulated heat can cause serious problems of overheating, which significantly degrades its efficiency, shifts the emission wavelength, and eventually reduces the lifetime.^{1} Thus, extracting the heat from LEDs attracts significant attention, especially using the materials with extremely high thermal conductivity.^{2,3} The cubic boron arsenide (BAs), a material with lattice thermal conductivity as high as ∼1500 W m^{−1 }K^{−1} at room temperature,^{4–7} has been synthesized, and its high thermal conductivity has been demonstrated recently. Its thermal conductivity is even comparable to diamond, the material with the highest thermal conductivity in nature.

With the advance in high thermal conductivity materials, the next bottleneck for heat dissipation is the thermal boundary conductance (TBC) between BAs and the LED wafer. A similar TBC issue also limits the applications of a diamond-based high thermal conductivity device.^{2,8} For practical applications in solid-state LEDs, the most commonly used low-cost gallium nitride based LEDs are grown on the substrate of silicon in the (111) plane.^{9,10} Thus, the TBC between the Si and BAs becomes extremely important to heat dissipation, the core of which is to understand the phonon-dominated thermal transport across the interface. While there has been significant progress in BAs single crystal synthesis very recently, there is no investigation on the TBC between the BAs and the other wafer materials, such as Si.^{11}

To predict the TBC across the BAs interface, molecular dynamics (MD) simulation is the most suitable method. This is because in high thermal conductivity materials, both the three-phonon scattering process and the four-phonon scattering process are important, although in most materials the four-phonon scattering process is negligible.^{12,13} This four-phonon scattering process is computationally very expensive in the first-principles related method. In contrast, within the MD simulations, this four-phonon scattering and other higher order scattering processes are automatically included,^{13,14} making them preferable methods for high thermal conductivity materials. It has been widely used to predict TBC across different interfaces. Generally, the nonequilibrium molecular dynamics (NEMD) method is preferable due to its simple configuration and direct measurement of the TBC from the temperature difference across the interface. Landry and McGaughey^{15} used the Stillinger–Weber (SW) potential and applied the NEMD method to predict TBC of the Si/Ge interface. They found that the TBC was about 320 MW/m^{2} K and was independent of the temperature below 500 K. The equilibrium molecular dynamics (EMD) is an alternative to consider. Gordiz and Henry^{16} conducted the EMD simulations and obtained the TBC in the range of 800 MW/m^{2} K.

In this study, we investigated the TBC between the BAs and Si across both (111) and (100) planes via the NEMD simulations where Si and BAs have a major phonon frequency mismatch. We also compared the TBC of the Si/BAs interface to the Si/Ge interface, where there was no significant phonon frequency mismatch and thus was deemed to have high TBC. This paper is organized as follows: In Sec. II, we describe the models and details for the density functional theory and MD calculations. In Sec. III, we first compare the accumulation thermal conductivity of cubic BAs, Si, and Ge as a function of frequency. Then, we test the validation of Tersoff potential used in this study before performing the MD simulations. Next, we calculate the length- and temperature-dependent TBC of the BAs/Si interface. At last, we explain the slight difference of TBC between BAs/Si and Si/Ge. Conclusions are drawn in Sec. IV.

## II. SIMULATION METHODS

### A. The first-principles calculation

We used the plane wave-based *Quantum Espresso* package^{17} and combined it with the Phono3py^{18} to obtain the force constants. Then, we used the Phono3py to calculate the thermal conductivity accumulation with respect to phonon frequency and the phonon relaxation time. We set the energy cutoff of 50 Ry, 3 × 3 × 3 supercell for force constants, and 31 × 31 × 31 k-mesh of wavevectors for thermal conductivity calculation. We first optimized the atomic structure. Then, we applied the finite displacement method to calculate the second-order harmonic and third-order anharmonic force constants. Then, with these force constants and based on the iterative solution to the linearized Boltzmann transport equation from phono3py package,^{18} we calculated the phonon scattering rates and thermal conductivity accumulation function with respect to phonon frequency.

### B. Molecular dynamics simulation

Both EMD and NEMD simulations are performed in this study. The former is used to calculate the thermal conductivity of the cubic BAs for the validation of the atomic interaction potential, while the latter is used to quantitatively extract the TBC in the BAs/Si interface for the comparison with that in the Si/Ge interface. In both simulations, all the MD were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package.^{19} Newton's equations of motion are solved by the velocity-Verlet integral algorithm with time step of 0.5 fs. We used the Tersoff potentials for Si and Ge.^{20} We adapted the required BAs atomic interaction potential from Benkabou *et al*.,^{21} which is also a convenient Tersoff potential. We tested this potential on its elastic, vibrational, and thermal transport properties in Sec. III B. The atomic interactions for the interface are obtained by a hybrid method.^{20}

In the EMD simulation, we first construct the lattice model of the cubic BAs with specified size, as shown in the inset of Fig. 1(a). Periodic boundary conditions are applied along three Cartesian directions. For each simulation case, we obtain the equilibrium structure using constant particle number, pressure, and temperature (NPT) simulations for 1 ns. Following the equilibration, we continue the EMD simulation under microcanonical ensemble (NVE) for another 5 ns. During this stage, the heat current of the simulation system is collected every 2.5 fs with the time step of 0.5 fs. The lattice thermal conductivity in the *a* direction (*a* = *x*, *y,* or *z*) is computed by the Green–Kubo formula,^{22–24}

where *V* is the volume of the simulation system, *T* is the absolute temperature, *k _{B}* is the Boltzmann constant,

*τ*is the correlation time, and

_{m}*J*

_{α}is the heat current component. $\u27e8J\alpha (\tau )J\alpha (0)\u27e9$ is the heat current autocorrelation function (HCACF) and

*τ*is the autocorrelation time. The

*τ*-dependent normalized HCACF for the cubic BAs at 700 K is shown in Fig. 1(a). The thermal conductivity is the integral with respect to

*τ*using Eq. (1), as shown in Fig. 1(b). In this work, we run three independent calculations to get the average thermal conductivity. Each run includes the average of all three directions. Thus, there are a total of nine calculations for the curve in Fig. 1. When the integral converges, its value is the lattice thermal conductivity.

In the NEMD simulation, we first constructed the atomic model of the BAs/Si interface as shown in the inset of Fig. 2. The heat baths (heat source and heat sink) are at the two ends. We selected the similar thicknesses for Si and BAs and varied the total thickness to check the length-dependent TBC. The periodic boundary conditions were applied parallel to the interface. After constructing the atomic model, we set the initial system temperature to 300 K in the MD simulation. Then, we applied a Berendsen heat bath of 300 K to all atoms for 0.5 ps to equilibrate the system. After that, we closed the thermal reservoir and ran the simulation for another 0.5 ns. Next, we fixed the temperatures of the heat source and heat sink to 330 and 270 K, respectively. After running the system for 1.5 ns, we obtained a steady temperature profile perpendicular to the interface as shown in Fig. 2. We obtained the temperature difference Δ*T _{inter}* across the interface by fitting the linear temperature regime inside both Si and BAs. The heat power across the entire cross section was obtained by fitting the accumulated adding power at the heat source or removing power at the heat sink. We averaged these two powers and referred it as $Q\u02d9$. With both $Q\u02d9and\Delta Tinter$, we can obtain the TBC

*G*by

_{inter}^{25}

where *A* is the cross-section area of the simulation model. We also evaluated the error of the calculated *G _{inter}* by exchanging the positions of the heat source and heat sink and running independent simulations with different initial conditions. We found that the uncertainty of

*G*is about 10%–20%. We also verified that using free or fixed boundary conditions along the heat flux direction would lead to similar

_{inter}*G*.

_{inter}## III. RESULTS AND DISCUSSION

### A. Thermal conductivity accumulation of BAs and Si

The thermal transport across the interface was affected by the overlap of phonons on both sides of the interface, especially phonon frequency and the phonon density of states (DOS).^{26–28} Thus, understanding the phonons’ contribution to thermal conductivity at different frequencies will be important to understand the TBC. As shown in Fig. 3, using the results from the above first-principles calculation, we obtained the phonons’ contribution to thermal conductivity at different frequencies for Si, BAs, and Ge. We first verified the results by calculating the thermal conductivity accumulation function with respect to phonon mean free path and validated it with the literature for Si.^{29} Our calculation only included the three-phonon Umklapp scattering and did not include the four-phonon scattering. Our accumulation function for BAs in Fig. 3 also agrees well with Lindsay *et al*.^{30} when we numerically integrate the thermal conductivity per frequency from 0 to 10 THz.

From Fig. 3, we found that the frequency range of major heat carriers is much narrower in BAs than that in Si. They do not have a major frequency overlap. For phonons with frequencies between 0 and 5 THz, they contribute to ∼70% of the thermal conductivity in Si, while only ∼20% in BAs, as shown in the light orange regime in Fig. 3. There is a significant frequency mismatch between Si and BAs normalized thermal conductivity accumulation function. This frequency mismatch can also be seen from the perspective of the phonon mean free path in the literature.^{4,29} Due to this frequency mismatch, the phonons in Si need other channels to talk to the phonons in BAs, such as high order phonon scattering^{31} or inelastic scattering.^{32} Thus, whether this frequency mismatch would lead to a low TBC would be an interesting and important topic. To further understand this frequency mismatch effect, we computed and compared the TBC of the Si/BAs interface to that of the Si/Ge interface in Sec. III C, in which the later interface has much less frequency mismatch as shown in Fig. 3.

### B. Validation of the Tersoff potential for BAs

The potential used in MD is extremely important because it plays the role in determining the results of the simulations. Thus, we tested the validation of the Tersoff potential for BAs before performing the NEMD to calculate the TBC. We first used the molecular simulation to validate the elastic properties, such as lattice constant *a*, the bulk modulus *B*_{0}, and the pressure derivative of the bulk modulus *B*′. We varied the lattice constant and calculated the size-dependent internal energy. Then, we fitted the internal energy by the Birch–Murnaghan state equation^{21,33} to find the lattice constant at its minimum, as shown in Fig. 4. We obtained the lattice constant of 4.777 Å, which is the same as in Ref. 21. We also obtained the bulk modulus *B*_{0} = 148 GPa and the pressure derivative of the bulk modulus *B′* = 4.7, which are close to the literature values *B*_{0 }= 141.1 GPa and $B0\u2032=4.1$.

We also calculated the phonon dispersion of bulk Si and BAs based on the Tersoff potential using our home-made codes, as shown in Fig. 5. The calculated dispersion for Si agrees well with the results of Lindsay *et al.*^{30} from first-principles calculation. It also works well with the acoustic branch for BAs. However, we also noticed that the optical mode for BAs is higher than first-principles calculation in the literature.^{30} This may be caused by the atomic interaction potential of BAs. For the optical mode, the adjacent atoms vibrate against each other. Thus, the interaction between the adjacent atoms is important for this mode. However, in the treatment of mixing rule,^{20,34} the parameters for B and As are averaged. Since the atomic mass difference between B and As atoms is very large (about 7), comparing to Si, this may introduce noticeable uncertainty for optical frequencies. Although there is a large uncertainty for optical modes, from the point of view of optical phonon group velocities and optical phonon-involved three-phonon scattering, these modes do not contribute to the thermal transport significantly. Regarding the group velocities, the optical phonons have low phonon group velocities as compared with that of acoustic phonons. Regarding the optical phonon-involved three-phonon scattering, there are two types of scattering processes, the aoo and aao processes. The aoo process is severely restricted in most materials due to the combined conservation conditions as showed by Lindsay *et al*.^{30} The aao process normally plays some role in scattering. However, the phonon bandgap between acoustic and optical branch in BAs is very large, which also prohibits the acoustic–optical scattering. Thus, the roles of optical modes are still small in the three-phonon scattering. However, it is admitted that the uncertainty of the four-phonon scattering caused by the optical phonon deviation cannot be excluded, especially at high temperature.^{12,13} In the future, better potentials need to be developed.

We further checked the BAs potential on its thermal transport properties. We calculated the temperature-dependent thermal conductivity using the EMD method as in Fig. 6. At all temperatures, our MD simulations slightly overpredicted the bulk thermal conductivity comparing to first-principles calculations.^{12} For example, the thermal conductivity at 700 K from the EMD simulation is about 402 W/m K, which is close to the prediction from the first-principles calculation.^{12} We also calculated the size-dependent thermal conductivities in the inset of Fig. 6, which show the cell size of 7 × 7 × 7, and periodic boundary conditions are sufficient to eliminate the size effect. One possible reason for the higher thermal conductivity from the present EMD is the missing of isotope scattering in our calculation, which may not be negligible even at high temperature in high thermal conductivity BAs.^{30}

### C. The TBC in BAs/Si interface

After validating the Tersoff potential for BAs above, we calculated the TBC for various interfaces as in Fig. 7 using NEMD simulations. We first built the atomic model for the Si(111)/BAs(111) interface as shown in the inset of Fig. 2. Using the same method, we also constructed atomic models for several other interfaces, such as Si(100)/BAs(100), Si(111)/Ge(111), and Si(100)/Ge(100). All the corresponding cross-section areas are 54.0 × 46.8 Å^{2}, 38.2 × 38.2 Å^{2}, 47.2 × 54.4 Å^{2}, and 44.5 × 44.5 Å^{2}. They are large enough to reduce the cross-section size effect when we apply the periodic boundary condition along the transverse directions for all MD models. The results are also comparable to or larger than the values of Si/Ge interfaces in the literature.^{35–37}

Regarding the NEMD method for TBC, we first need to consider the simulation length effect in the heat flux direction. As shown in Fig. 7, we did the simulation at three or four lengths, close to 15 nm, 30 nm, and 60 nm. We did not find a significant TBC enhancement when the length was doubled or quadrupled from 15 nm. This means the length 30 nm is sufficient to eliminate the length effect. This length dependence agrees with the literature^{16} on Si/Ge. We obtained the TBC for the Si/BAs (100) plane as ∼300 MW/m^{2} K. It is a relatively high TBC value, as we compare the Si/BAs interface with other interfaces and consider there is a significant frequency mismatch between Si and BAs. This seems to contradict with the thermal conductivity accumulation function in Fig. 3.

To answer the question of how the phonon frequencies mismatch in both sides of the interface affects the TBC, we also calculated the TBC across the Si/Ge interface along both (111) and (100) planes. From Fig. 3, we can see that the normalized accumulation function of Si and Ge is very close to each other. This means a large number of phonons in Si can transport to Ge at the same frequency. Thus, the TBC in the Si/Ge interface should be high. However, when we compare the TBCs of Si/BAs and Si/Ge on both (111) and (100) planes as shown in Fig. 7, we found there is no significant difference between the TBCs on Si/BAs and Si/Ge in the (100) or (111) plane. The TBC of the Si/BAs interface is only ∼40% lower than the Si/Ge interface in the (111) direction. All the TBCs are similar ranging from 200 to 400 MW/m^{2} K. This means the frequency mismatch on both sides does not significantly reduce the TBC. It also indicates that the Si/BAs interface is efficient for heat dissipation, which may have important applications in thermal management.

We also changed the simulation temperature and obtained the temperature-dependent TBC of the BAs/Si interface. Since there only exists slight difference between the TBCs of BAs(111)/Si(111) and BAs(100)/Si(100) interfaces, we only present the TBC of the Si(111)/BAs(111) model to save the computational cost. Figure 8 shows that the TBC presents a slight increase when the temperature increases from 300 to 700 K. This is reasonable because the estimated Debye temperature of BAs is above 580 K from the cut-off acoustic phonon frequency. Most acoustic phonons in BAs have been excited at room temperature. Therefore, the TBC of BAs/Si is not sensitive to temperature above room temperature.

### D. High TBC across BAs/Si interface

To explain the high TBC value across the Si/BAs interface, we calculated the phonon density of states (DOS) for Si and BAs as shown in Fig. 9(a). The DOS is obtained from the Fourier transform of the velocity autocorrelation function with additional MD simulation.^{25} We can see in the light orange regime of 5–8 THz, which is the dominated phonon frequency for heat conduction in BAs, that the DOS has a big overlap between Si and BAs for the acoustic phonons. This DOS overlap indicates the possible high probability of thermal transport across the interface.^{38}

To explain the frequency mismatch in Fig. 3, we also calculated the phonon scattering rate from first-principles, as shown in Fig. 9(b). We can see that the major low acoustic phonons scattering rates are in the frequencies between 5 and 8 THz. These low scattering rates explain why the major heat carriers in BAs have the frequency between 5 and 8 THz. However, these scattering rates of acoustic phonons in Si mostly fall in the frequency range lower than 5 THz.^{29} This explains the frequency mismatch at the beginning. Thus, it is the *mismatch of scattering rates that determines the frequency mismatch in thermal conductivity accumulation function* in Fig. 3; and it is the *DOS overlap that determines the high TBC* in Fig. 9.

We also compared the TBC of Si/BAs with other interfaces, such as Pb/Diamond, Al/Al_{2}O_{3}, or Al/TiN,^{39} as shown in Fig. 10. The TBC value is relatively high compared with these interfaces. This means the BAs is an effective heat conductor as a single crystal itself and also as an interface material with Si. There may be other possible reasons for this high TBC, such as new interface mode.^{16} Gordiz and Henry^{16} showed that the interface modes, which have high frequencies, also played important roles in the thermal transport across the Si/Ge interface. There are some other explanations of these interface modes, such as high order phonon scattering,^{31} inelastic scattering,^{32} etc. All these explain the high TBC across the interfaces. Thus, the high TBC in the Si/BAs interface indicates that BAs is an excellent material for heat dissipation and will have broad applications in thermal management in LEDs and other high-power electronics.

## IV. CONCLUSIONS

In summary, we predicted the thermal boundary conductance of the Si/BAs interface using the NEMD simulations. We found that the TBC is as high as 200–300 MW/m^{2} K even though there is a large phonon frequency mismatch of the thermal conductivity accumulation function between the Si and BAs. We also found the lattice directions, such as the (100) and (111) directions, do not affect the TBC. Comparing with the Si/Ge interface, we found that the TBCs are not significantly different even though there is a frequency mismatch in the Si/BAs interface. Through calculating the phonon DOS and scattering rate, we found that the frequency mismatch in the thermal conductivity accumulation function is caused by the frequency mismatch of the phonon scattering rate, and the overlap of acoustic phonon DOS causes the high TBC. Additionally, we expect that the interface mode may also play an important role across the interface.

## ACKNOWLEDGMENTS

Zhiyong Wei acknowledges financial support from the National Key R&D Program of China (No. 2017YFB0406000) and Southeast University “Zhongying Young Scholars” Project. Fan Yang acknowledges start-up support from the Stevens Institute of Technology. The authors also acknowledge help and support from Maria Chan and Subramanian Sankaranarayanan from the Center for Nanoscale Materials, an Office of Science user facility, and were supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences under Contract No. DE-AC02-06CH11357.