High throughput substrate screening for interfacial thermal management of β-Ga₂O₃ by deep convolutional neural network

In the rapid development of nano- and micro-electronic devices, thermal management is a crucial component to prevent overheating, which might lead to failure, wearout, and the lifetime shortening of the devices. Thermal transport across the interface is a critical part of the thermal dissipation process in those devices.^1–3 The performance of thermal transport between two materials in contact can be measured by the interfacial thermal conductance (ITC) through ITC = Q/ΔT, where Q and ΔT are heat flux and temperature difference between the two materials at the interface, respectively.^1–4 As miniaturization continues, in particular, in nanostructured devices, the characteristic length scales of devices are much shorter than the intrinsic mean free path of phonons in the counterpart bulk materials, and thus, ITC becomes the most important factor for thermal transport. In reality, thermal transport at the interface can be affected by several factors such as binding energy, roughness, and the presence of defects and impurities that could occur by the mixture of atoms at intermediate layers at the interface.⁵ Phonon reflections caused by the acoustic properties of the adjacent materials take place across the interface even if the two materials are in perfect contact. Therefore, various parameters control ITC, which makes it even more difficult to predict or approximate.^1–6 

For more accurate ITC computations, more computationally expensive calculations must be implemented such as nonequilibrium molecular dynamics (MD) simulations that naturally consider both harmonic and anharmonic phonon scatterings^12,13 and Green's function method where the phonon anharmonicity can even be considered.^14,15 Despite the shortcomings of DMM compared to experimental results,¹⁶ computational results by DMM are still useful to capture the tendencies in thermal transport across interfaces.¹⁷ More importantly, the DMM offers a fast-screening approach due to the less expensive computations before experimental synthesis and measurements. However, prediction using the DMM model requires detailed information about the mode level phonon properties such as phonon group velocity and phonon DOS, both of which are not easy to obtain in a high-throughput manner.

The advent of machine learning (ML) and artificial intelligence (AI) has revolutionized many aspects of modern science and technology and has sparked significant interest in the thermal science community in recent years. As a means to make computations even faster, ML algorithms are utilized on data generated by density functional theory (DFT) calculations since they can be utilized to predict various material properties such as mechanical,^18–20 thermal,^21–23 magnetic,²⁴ and optical²⁵ properties. In this work, we first performed high-throughput DFT calculations on comprehensive phonon properties of 1633 heat sink materials and then calculated their ITC when being contacted with β-Ga₂O₃ as a heat source by the DMM model. Here, the β-Ga₂O₃ was chosen as heat source as a case study since it has recently gained significant scientific interest for future power devices because of its inherent material properties such as extremely wide bandgap, outstanding Baliga's figure of merit, large critical electric field, etc. β-Ga₂O₃ has been widely used in high-electron mobility transistor (HEMT) applications.^26,27 Using DFT-ITC results as training data, we trained our recently developed convolutional neural network (CNN) model that utilizes the fused orbital field matrix (OFM) and composition Magpie descriptors and then deployed the trained CNN model to predict ITC of large-scale new heat sink materials in contact with β-Ga₂O₃. This work can be helpful in inverse designing new substrates for HEMT devices that utilize β-Ga₂O₃ as a heat source.^28,29

Figure 1 shows the schematic chart for the entire workflow conducted in this study to accelerate the search for promising substrates or heat sink materials for heat dissipation of β-Ga₂O₃. The first step is to obtain the structures from the open quantum material database (OQMD), and we re-optimized them with our own DFT parameters (see computational details of DFT calculations below). Then, the second- and third-order interatomic force constants (IFCs) of the majority of structures were obtained from our previous works^{21–23,30–32} with a small amount of additional DFT calculations utilizing the same compression sensing lattice dynamics (CSLD) method to calculate the phonon dispersions. The second-order IFCs confirm the dynamic stability through the absence of negative phonon frequencies in the Brillouin zone. Then, almaBTE calculations using both second- and third-order IFCs of materials as input are performed to create a dataset for β-Ga₂O₃/X interfaces having ITC as target (here “X” denotes a substrate or heat sink material in contact with β-Ga₂O₃). The dataset is then used to train our developed fused OFM + Magpie CNN model, and subsequently, the model is used to make ITC predictions on new candidate substrates in contact with β-Ga₂O₃. For the screening pool, we prepared 32 716 structures with non-zero bandgap as a special requirement for thermal management applications. All 32 716 structures were screened out from the OQMD database and have been successfully re-optimized by DFT calculations with our own computational parameters. The structures that were not successfully optimized were discarded.

All the 1633 crystal structures used for training our developed fused OFM + Magpie CNN model are downloaded from the OQMD database in crystallographic information file (CIF) format.³³ All the structures were re-optimized by first-principles DFT calculations utilizing Vienna Ab-Initio Simulation package.^34–36 For re-optimizing the structures, we use their primitive cells with total energy and force convergence criteria of 10⁻⁸ eV and 10⁻⁴ eV/Å, respectively. The atomic positions, cell shape, and sizes of the primitive cells are all allowed to change in the re-optimization process to reach the minimum energy in each structure. The Perdew–Burke–Ernzerhof (PBE)³⁷ is selected from the generalized gradient approximation (GGA) to characterize the exchange-correlation effects of electrons within the linear projector augmented wave (PAW) generalization of the pseudopotentials method³⁸ to allow DFT to perform calculations with higher efficiency. The kinetic energy cutoff of 520 eV was selected to limit the number of plane waves utilized as basis functions at each wavevector in the k-space. The Monkhorst–Pack k-mesh was applied to sample the Brillouin zone.³⁹ The number of k-points for electrons was determined such that the multiplication between the lattice vector and the number of k-points is at least 60 at each lattice vector to ensure high-quality k-mesh. After primitive cells were re-optimized, supercells were generated based on the optimized primitive cells, in which each atom was randomly displaced by a displacement of 0.03 Å. For each primitive cell, 16–30 supercells were generated, with a detailed number of supercells depending on the symmetry of the materials. The atomic forces of the supercells were then evaluated by VASP calculation with energy and force convergence criteria of 10⁻⁶ eV and 10⁻⁴ eV/Å, respectively. After that, the second-order (harmonic) and third-order (anharmonic) interatomic force constants (IFC) were generated using the compressive sensing lattice dynamics (CSLD) method,^40–42 and the IFCs were output in Phonopy⁴³ format. ShengBTE package was then utilized to obtain the lattice thermal conductivity (LTC)⁴⁴ through an iterative method. The LTC of β-Ga₂O₃ was calculated by our DFT + BTE to be 11.6, 20, and 16.6 W/mK in the [100], [010], and [001] direction, respectively, which are in good agreement with those reported in Ref. 45. With second- and third-order IFCs, the almaBTE package⁴⁶ was utilized to approximate the interfacial phonon transmission coefficient using the DMM model and subsequently obtain ITC through Monte–Carlo simulations. In the almaBTE calculations, two slabs composed of β-Ga₂O₃ as the heat source and another substrate as the heat sink were set up, with each slab of 100 nm thick. It should be noted that, since β-Ga₂O₃ is a highly anisotropy material, the ITC would depend on the specific axis or orientation chosen when running almaBTE. After testing all three cartesian axes, we find that the heat transport along [010] cartesian axis yields the highest ITC, which can be understood in terms of the highest LTC in the same direction. Therefore, all ITC results presented in this work refer to the orientation of [010] for β-Ga₂O₃. For the phonon wavevector q in the almaBTE calculations, 50 q-point line density for each lattice vector was implemented. The simulation temperature difference between the heat source and heat sink is 5 K. The number of particles in the simulation is 1 × 10⁶, which is large enough after testing. The temperature is output on 300 equally spaced bins sliced perpendicularly to the heat flux direction in the two slabs. The ITC between β-Ga₂O₃ and the substrate is calculated by the temperature difference at the interface.

Figure 2 shows the schematic chart for our recently developed fused orbital field matrix (OFM) + Magpie CNN model. OFM⁴⁷ is based on the representation of the valence shell electrons of a local chemical environment by considering the sum of the weighted vector of all atoms in that environment. More details on the OFM calculation method can be learned from Ref. 47. Magpie features⁴⁸ are statistical results of maximum, minimum, mean range, and mode of the elemental descriptors generated based on the composition of the material. OFM and magpie descriptors do not require expensive computations to acquire them. Matminer⁴⁹ can generate such descriptors for thousands of materials in less than 1 s. The layer numbers shown in red color are for both OFM⁴⁷ and Magpie⁴⁸ descriptors. The OFM we used is a 32 × 32 matrix or second-rank tensor. Such matrix can be used as an input for a 2D CNN as seen in the conv2D in Fig. 2 layer 1 and layer 2(b) shown in the texts in red color. The feature matrix from the two separate path conv2D in layer 2(a) and layer 2(b) is combined and proceeds to go to conv2D in layer 3 and then max pooled before going through another conv2D shown indicated by the red number 4. The number of kernels in OFM forward propagation conv2d layers is as follows from the red numbers in the “ResidualLikeBlock”: 1–5, 2a–3, 2b–3, 3, 4–3. The number of filters in OFM forward propagation Conv2d layers is as follows: 1–32, 2a–32, 2b–32, 3–64, 4–64. In the magpie forward propagation Conv2D layers, the number of kernels is as follows: 1–3, 2a–3, 2b–3, 3–3, 4–3, and the number of filters is as follows: 1–32, 2a–48, 2b–48, 3–64, 4–64. The ReLU activation function is used in all the conv2d layers and dense layers except for the last dense layer that has no activation function. The padding is the “same” in all the conv2D layers. The pool size is (2, 2) in the MaxPooling2D layer. The 2D feature matrix is then flattened to become a feature vector that has the shape of 16 384. The Magpie descriptors can be represented as a 132-feature vector and can be shaped into a 12 × 11 matrix to be used as an input to a conv2D. The Magpie matrix follows a similar path of layers before flattening and becomes a feature vector that has the shape of 2304. The flattened feature vector from the OFM path and flattened feature vector from the Magpie path are then concatenated before progressing to three dense layers to make a prediction on ITC of substrates in contact with β-Ga₂O₃. The model is to train the 1633 materials dataset. The 1633 materials dataset is split into 80%/10%/10% for training, validation, and testing sets, respectively. The learning rate used in the AdamW optimizer is 0.001. A similar model was used in this work.⁵⁰ However, the model was used for classification not regression as in this work. The entire model flowchart is built on Tensorflow.⁵¹ 

Figure 3(a) shows the training and validation losses with respect to the number of training epochs. The loss function in each epoch is the mean absolute error (MAE) loss function as implemented from Tensorflow. The maximum number of epochs is 60 epochs, which is more than enough to reach the lowest loss. As a matter of fact, the model plateaus do not significantly reduce the loss after the tenth epoch, which means the model training converges very fast. In Fig. 3(b), we compare the ITC in base 10 logarithm values between full DFT calculations and prediction by the fused OFM + Magpie CNN model with R² score labeled. Generally speaking, the R² score of an ML model prediction should be as high as possible (e.g., close to unity) since the metric denotes the consistency of the results from testing or the model’s previously unseen dataset of ITC values compared to the predictions by the model. The R² score of our trained fused OFM + Magpie CNN model is as high as 0.874, which confirms that our fused OFM + Magpie CNN model is well trained with high quality and can, thus, be utilized to predict ITCs for new or previously unseen substrates in contact with β-Ga₂O₃.

Since there are many previously widely used ML models available, then it is worth comparing our fused OFM + Magpie CNN model with those models to affirm its superiority in ITC predictions. To this end, our regression results from the fused OFM + Magpie CNN model are compared with two well-known and traditional ML models: (1) gradient boosting (GB)⁵² and (2) random forest (RF).⁵³ Both ML models used the Magpie descriptors to train their respective models. These two ML algorithms are imported from sci-kit learn library.⁵⁴ The comparison among the three models is implemented by comparing the metrics among the models, i.e., R² score and MAE. Table I shows the comparison results of the metrics among the three trained models. The R² score of our developed fused OFM + Magpie CNN model is the highest among three ML models. MAE should be as low as possible to substantiate that model's prediction (i.e., predicted ITC) is not far off from the true value (i.e., “true” ITC). We find that the MAE of our developed CNN model is the lowest among the three models. From these two metrics, we can conclude and affirm that our developed CNN model has superior performance and tremendous capability compared to the other two traditional ML algorithms in predicting ITC of substrates in contact with β-Ga₂O₃.

With the well-trained fused OFM + Magpie CNN model, we continue to deploy the model to predict the ITC of 32 716 structures in contact with β-Ga₂O₃. Figure 4 shows the statistical distribution of the ML dataset of 1633 DFT materials and 32 716 materials in the prediction dataset in terms of the number of atoms, mass density, and volume of primitive cells. The prediction pool is much larger and contains a wider range of structures compared to the ML dataset that was used to train, validate, and test the ML model. Here, we need to emphasize that the ITCs in base 10 logarithm values were first predicted, and then they were converted back to the normal values in units of MW/m²K. The top 10 structures with the highest predicted ITCs as promising substrates for cooling β-Ga₂O₃ are presented in Tables II and III, respectively. Relevant structure information such as formula, bandgap, formation energy, etc., are also provided. Interestingly, all top 10 cubic substrates happen to have the same space group number of 227 and the same prototype AB₂C₄, although our screening pool includes lots of materials with other space group numbers such as 225, 221, and 216, and other prototypes as well such as AB, ABC, ABC₂, ABCD, etc. This means similar atomic structure and other similar structural features will likely lead to similar ITC, as we will see shortly from the analysis of DFT-ITC data. In contrast, the top 10 noncubic substrates have slightly diverse space groups and prototypes, primarily because there are many more space groups and prototypes in the screening pool. Again, similar material families in particular with the same prototype and close average atomic mass tend to have similar ITCs.

Since we have prepared many Magpie atomic descriptors for all structures in our training dataset during traditional ML training, it is then intuitive to calculate the Pearson correlation coefficients⁵⁵ of ITC with those Magpie atomic descriptors. Pearson correlation values range from −1 to +1. Positive (negative) values denote a directly proportional (inverse) correlation between two variables. The high value of Pearson correlation regardless of the sign (i.e., close to +1 or −1) indicates a strong correlation between the two variables, while the Pearson correlation coefficient close to zero, regardless of the sign, denotes a weak or no correlation between two parameters. The three Magpie descriptors selected to show in Table IV are the mean atomic number, mean atomic weight, and mean ground state (GS) volume per atom (PA), which have the highest Pearson correlation coefficients of −0.4995, −0.4987, and −0.4637, respectively, among the total 132 Magpie descriptors we have studied. The coefficients are represented in absolute value because the strength of the correlation coefficient is one of the crucial factors to correlate with ITC regardless of the sign. The reason why these descriptors have high Pearson correlation coefficients is because the constituent species in β-Ga₂O₃ all have medium to low atomic weights, which means the material possesses phonon states at high frequencies in its phonon dispersion in the reciprocal space.^56,57 As implied in Eqs. (1) and (2), similar vibrational frequencies between the heat source and substrate are one of the crucial factors in acquiring high ITC. Therefore, materials with medium to low atomic weight or atomic number should have high frequencies and possibly high phonon DOS overlap with β-Ga₂O₃. Regarding the high Pearson correlation coefficient for the mean GS volume per atom, it is due to the light elements having low GS volume per atom, and those light elements also tend to have high frequencies in their phonon dispersions. Therefore, the mean GS volume per atom descriptor happens to have the same trend in the periodic table as the atomic weight, which made its correlation high with ITC. The same analysis for the mean GS volume per atom can be applied to other descriptors such as covalent radius because it follows the same periodic table trend.

From Table IV and Fig. 5, it is confirmed that ITC is inversely proportional to the similarity in atomic number or atomic mass, which comes from the vibrational similarity or phonon DOS overlap between β-Ga₂O₃ and the substrate. Therefore, we continue to quantitatively investigate the relationship between ITC and phonon DOS overlap. Figure 6 manifests ITC values of β-Ga₂O₃ with other substrates along with phonon DOS between β-Ga₂O₃ and those substrates. LTC of the substrates is added as a color map to the analysis since it adds a measure for the phonon group velocity of the substrate, which is also another crucial parameter for ITC. The clear and obvious trend that can be seen from Fig. 6 is that ITC is equal to zero if phonon DOS overlap between β-Ga₂O₃ and another substrate is also equal to zero, which is well expected under the assumption of the DMM model used by the almaBTE package. This analysis and results are completely understandable since Eqs. (1) and (2) reveal that zero phonon DOS overlap should have zero ITC as well. In other words, phonon DOS overlap must not be zero in order to make a nonzero or finite ITC. It is also shown from the bottom-right region in Fig. 6 that, although phonon DOS overlap is a bit high between 0.3 and 0.6, the ITC is relatively low. This is because the LTCs of most substrates are low (the corresponding LTC color is dark blue). That could also convey that those substrates have low group velocity and as a result low ITC. In that same range (i.e., phonon DOS overlap from 0.3 to 0.6), the data higher than the bottom-right (i.e., higher ITC) have higher LTC as indicated by the color in the LTC color bar. This indicates that higher LTC leads to higher ITC for the same phonon DOS overlap, which can be explained by the fact that the higher LTC usually has high group velocity and consequently higher ITC. However, this trend is not always true. In Fig. 6, we find that there are lots of substrate materials with very high ITC when being contacted with β-Ga₂O₃, but those materials do not necessarily have high LTC. Specifically, for ITC above 800 MW/m²K, many substrates only have moderate to low LTC values, roughly in the range from 10 to 30 W/mK. The high ITC is then primarily resulted from the high phonon DOS overlap. Overall, the results in Fig. 6 demonstrate that for heat sources with low LTC like the case of β-Ga₂O₃ (the LTCs in three directions are between 10 and 20 W/mK obtained by our DFT calculations), the phonon DOS overlap between two neighboring materials is much more important than the LTCs of substrates themselves in terms of achieving high ITC at the interfaces.

Since the high phonon DOS overlap is a necessary condition to obtain high ITC between a heat source and substrate as understood from Eqs. (1) and (2) and visually portrayed in Fig. 6, particularly for two contacting materials with low LTC, in Fig. 7, we show two representative examples, (a) namely, MgO and (b) YN, of the phonon DOS overlap with β-Ga₂O₃. These two substrates were chosen because their ITCs are the highest among all 1633 DFT-ITC datasets. Specifically, the ITC for β-Ga₂O₃/MgO and β-Ga₂O₃/YN interface is as high as 1211 and 1224 MW/m²K, respectively. It can be seen from Fig. 7 that high phonon DOS overlap exists between β-Ga₂O₃ and the two substrates, with the calculated values of 0.56 and 0.63 for β-Ga₂O₃/MgO and β-Ga₂O₃/YN interface, respectively. The phonon DOS overlap spans acoustic and optical phonon branches. Moreover, the phonon group velocity seems to be high in MgO in the acoustic branch and most of the phonon modes in the optical branch except for those from 15 to 17.5 THz, as shown in Fig. 7(c). A similar phenomenon is observed in YN where high phonon group velocity is observed in Fig. 7(d) for all phonon modes, particularly for some optical phonon modes with frequency above 7.5 THz.

Although Fig. 7 shows the phonon modes that possibly contribute to transporting heat from the phonon DOS overlap and the high group velocity phonon region that might transport more heat across the interface, it does not quantitatively tell us which phonon modes can actually transport how much heat across the interface, since transmission coefficient is another very important factor to interfacial thermal transport according to Eqs. (1) and (2). Besides, almaBTE considers several scattering mechanisms such as intrinsic phonons scattering, interface scattering, and absorption scattering. Therefore, even with phonon DOS overlap and high group velocity, the scattering mechanisms implemented in almaBTE might hinder transporting higher heat fluxes. To this end, Fig. 8 shows which phonon modes and frequencies transport how much heat. For the β-Ga₂O₃/MgO interface shown in Fig. 8(a), phonon frequencies from around 0.5 THz to roughly 8 THz transport the highest amount of heat compared to other phonon modes that overlap between 8 and 17 THz. That can be explained by the high group velocity in the acoustic branches and less phonon scattering compared to the higher frequency phonon modes in the optical branches, which might have high phonon group velocity in most of the states but high phonon scatterings as well. For β-Ga₂O₃/MgO interface shown in Fig. 8(b), high heat flux is transported by the phonons with frequency between 0 and 9 THz, which indicates the higher group velocities at those frequencies and less phonon scattering. Similar to the β-Ga₂O₃/MgO interface, although there is a large overlap between those two materials for phonon states from 9 to 17 THz, these phonon modes do not contribute much to transporting a significant amount of heat flux, which is mainly due to high phonon scattering at those phonon states. Overall, the spectral heat flux analysis in Fig. 8 demonstrates that the high phonon DOS overlap and high phonon group velocity along with low scattering rates are required for achieving high ITC.

In summary, the objective of this work is to accelerate screening suitable substrates for cooling β-Ga₂O₃, a common material used in HEMT devices. The potential substrate candidates should have high ITC in contact with β-Ga₂O₃ to dissipate as much heat as possible to prevent the devices from wearing out and to prolong their lifetime. To this end, the ITC of 1633 substrates in contact with β-Ga₂O₃ is calculated using almaBTE with the required second- and third-order IFCs computed directly from high-accuracy DFT calculations. DFT-ITC results reveal that the phonon DOS overlap plays an important effect on the ITC, and the phonon DOS overlap of about 0.3 or higher is a necessary condition to achieve high ITC. Group velocity plays another crucial role in attaining high ITC. For phonons in the finite (non-zero) DOS overlap region, the phonon group velocity becomes one of the dominant factors for ITC. The third factor affecting ITC is the scattering at the interface. Even for the phonon modes at various frequencies that overlap between β-Ga₂O₃ and a specific substrate, it is not guaranteed that high heat flux can be transported across the interface by those modes since they might have strong scatterings. The β-Ga₂O₃/YN and β-Ga₂O₃/MgO interfaces possess the highest ITC of 1224 and 1211 MW/m²K, respectively, among the 1633 DFT-ITC datasets. Both substrates fulfill the three governing factors or requirements mentioned above, namely, high phonon DOS overlap, high group velocity, and weak scatterings. With 1633 DFT-ITC datasets, we trained our recently developed convolutional neural network model that utilizes the fused orbital field matrix and composition descriptors. Our model proved to be superior in performance to traditional machine learning algorithms such as random forest and gradient boosting. We further deployed the CNN model to predict the ITC of 32 716 structures with non-zero bandgap in contact with β-Ga₂O₃. The CNN model predicted the top 20 cubic and noncubic substrates with ITC on the same level (around 1000 MW/m²K) as Ga₂O₃/YN and β-Ga₂O₃/MgO interface. We also performed Pearson correlation analysis and found three Magpie descriptors that have a strong correlation with ITC, namely, mean atomic number, mean atomic weight, and mean ground state volume per atom. These descriptors based on the primitive cells of substrates are easy to calculate as soon as the structures are optimized, which provide a new route for quickly screening potential substrates from large-scale material pools for high-performance interfacial thermal management of HEMT devices.

Research reported in this publication was supported in part by NSF (Award Nos. 2110033 and 2311202), SC EPSCoR Program under Award No. 23-GC01, and an ASPIRE grant from the Office of the Vice President for Research at the University of South Carolina (Project No. 80005046).

The authors have no conflicts to disclose.

Mohammed Al-Fahdi: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Ming Hu: Conceptualization (lead); Formal analysis (supporting); Funding acquisition (lead); Project administration (lead); Resources (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. The fused OFM+magpie CNN model code is available free of charge at https://github.com/Mofahdi/OFM-magpie-CNN-Model-for_ITC_prediction. The almaBTE scripts and tutorials are available free of charge at https://github.com/Mofahdi/almaBTE_io.