Dissimilar thermal transport properties in κ-Ga₂O₃ and β-Ga₂O₃ revealed by homogeneous nonequilibrium molecular dynamics simulations using machine-learned potentials Open Access

Ultra-wide bandgap semiconductors, such as $Ga 2 O 3$⁠, diamond, and AlN, have also become the focus of attention materials for next-generation electronics and optoelectronics. Owing to a bandgap of about 5 eV, exceptional breakdown electrical field, and cost-effective production, $Ga 2 O 3$ offers considerable promise for ultrahigh power device applications.^1–3 However, the lattice thermal conductivity (LTC) of $Ga 2 O 3$ is subpar, leading to pronounced heat dissipation issues in certain semiconductor devices.¹ Understanding the phonon thermal transport in $Ga 2 O 3$ is crucial for its practical applications.

The $Ga 2 O 3$ crystal actually exists in five distinct phases: $α$⁠, $β$⁠, $γ$⁠, $σ$⁠, and $ε$ (sometimes referred to as $κ$⁠). Of these, the $β$ phase (space group C2/m) is the most stable one, which has been extensively explored for applications in deep-ultraviolet transparent conductive electrodes,⁴ solar blind detectors,^5,6 high-performance field effect transistors,⁷ Schottky rectifiers,⁸ and high temperature gas sensors.¹ In recent years, further efforts have been made to overcome the poor thermal stability and immature synthesis methods of other phases. Hexagonal crystal $ε− Ga 2 O 3$ (space group P63mc) was reported to be the second most stable phase obtained at low temperatures,^9,10 whereas subsequent studies claimed that the crystal structure of the polycrystalline form at low temperatures has an orthorhombic structure at the nanoscale (5–10 nm), named $κ− Ga 2 O 3$ (space group Pna21).^11,12 The strong polarization in $κ− Ga 2 O 3$ is a prominent feature that $β− Ga 2 O 3$ does not possess, which may benefit potential device applications. For example, $κ− Ga 2 O 3$ can form heterojunctions with other semiconductors such as GaN or AlN, and its polarization was utilized to regulate interface transport while effectively alleviating thermal problems.^13,14 Considering the distinct crystal structures of $β$ and $κ− Ga 2 O 3$⁠, a comparative understanding of the LTC for both phases is essential for their device thermal design, particularly as the $κ$ phase remains unexplored.

Currently, growing large-scale single-crystalline $κ− Ga 2 O 3$ is a challenging endeavor experimentally, and assessing its LTC is anticipated to be even more demanding. From a computational perspective, its LTC can be predicted using atomistic simulation techniques, such as molecular dynamics (MD) simulations^15–17 and the combination of Boltzmann transport equation (BTE) with the anharmonic lattice dynamics (ALD) method.^18–20 While atomic interactions can be derived from either empirical potentials or quantum mechanical density functional theory (DFT) calculations, intricate crystals like $κ− Ga 2 O 3$⁠, which possess 40 atoms in their primitive cell, introduce substantial hurdles—accuracy concerns for empirical potentials and computational overhead for DFT. Recently, the prediction of thermal transport properties in materials has been significantly accelerated through the combined use of machine learning approaches and atomistic simulation methods.²¹ Specifically, machine-learned potential (MLP)-based large-scale MD simulations have been demonstrated to be a reliable approach to calculate the LTC of complex crystals including amorphous silicon,²² amorphous $HfO 2$⁠,²³ amorphous silica,²⁴ metal-organic frameworks,²⁵ and violet phosphorene,²⁶ which can account for phonon anharmonicity to arbitrary order. Several MLPs for $Ga 2 O 3$ have been developed, addressing both the perfect bulk system^15,19,27 and more intricate, disordered structures.^28,29

In this work, we apply the neuroevolution potential (NEP) framework^30–32 to develop two MLPs on demand for $Ga 2 O 3$ against quantum-mechanical DFT calculations, one for $β$ phase and one for $κ$ phase. We choose the NEP approach here, because this method has been demonstrated to be highly efficient.³² We apply the developed NEPs to perform extensive homogeneous non-equilibrium molecular dynamics (HNEMD) simulations to investigate the LTC of the two phases of $Ga 2 O 3$⁠. Our results show that the LTC of the $κ$ phase is much lower than that of the $β$ phase and exhibits low anisotropy and weak temperature dependence.

All DFT calculations were performed using projected augmented wave method³³ with Perdew–Burke–Ernzerhof functional of generalized gradient approximation³⁴ implemented in the VASP package.^35,36 The kinetic energy cutoff was set to 520 eV, and the convergence value for the total energy was $1× 10 − 6$ eV. The stopping criteria for structural optimizations were that the maximum residual Hellmann–Feynman force on atoms was less than $1× 10 − 3 eV/Å$⁠. A $10×10×6$ and $10×6×6$ $Γ$-centered k-point grid was employed for the primitive cell of $β− Ga 2 O 3$ and $κ− Ga 2 O 3$⁠, respectively. Additionally, to validate the accuracy of the MLP, the phonon dispersions of $β− Ga 2 O 3$ and $κ− Ga 2 O 3$ were further calculated on $2×2×2$ and $3×1×1$ supercells, respectively, using the density functional perturbation theory together with the Phonopy package.³⁷ 

Meanwhile, according to widely used conventions, we have optimized the corresponding conventional unit cell as the initial cell for the MD simulation. The $κ− Ga 2 O 3$ lattice structure remains unchanged, but the number of atoms in the conventional unit cell of $β− Ga 2 O 3$ is higher than that in the primitive cell. The calculated lattice constants of $β− Ga 2 O 3$ are $a=12.468$⁠, $b=3.087$⁠, and $c=5.716$ Å, while these values for $κ− Ga 2 O 3$ are $a=5.074$⁠, $b=8.703$⁠, and $c=9.309$ Å, which are close to the reported values.^12,38 As shown in Fig. 1, the conventional unit cell of $β− Ga 2 O 3$ has 20 atoms, while $κ− Ga 2 O 3$ has a relatively large number of 40 atoms. In the calculations of the reference datasets for NEP training, the $2×2×2$ and $2×1×1$ supercells were used for $β− Ga 2 O 3$ and $κ− Ga 2 O 3$⁠, respectively, with $1×4×2$ and $3×3×3$ k-point grids.

We used the identical method to construct the reference datasets for the $β$ and $κ$ phases of $Ga 2 O 3$⁠. The reference structures were obtained by ab initio molecular dynamics (AIMD) simulations and random perturbations. For each phase, the AIMD simulations were run under an isothermal ensemble with temperature increasing linearly from 10 to 1000 K for 10 000 steps and a time step of 1 fs. We uniformly extracted 1000 structures from AIMD simulations for each case. For random perturbations, 200 structures were generated with random cell deformations ranging from $−$4% to 4% and atomic displacements less than 0.1 Å, all based on optimized structures. For both phases, our total dataset has 1200 structures including energy, atomic forces, and virials obtained from DFT calculations as outlined in Sec. II A. We randomly divided the total dataset into a training set of 1000 structures and a test set of 200 structures.

After obtaining the training set and test set, we applied the third generation of the NEP framework³² implemented in GPUMD (version 3.5) to train the MLP for $β$ and $κ$ phases of $Ga 2 O 3$⁠, which were denoted as NEP- $β$ and NEP- $κ$⁠, respectively. NEP used a feedforward neural network to represent the site energy with atomic cluster expansion (ACE)³⁹-like descriptor components including radial and angular terms. The parameters of NEP models were optimized using the separable natural evolution strategy (SNES),⁴⁰ with the loss function defined as a weighted sum over the root mean square error (RMSE) values of energy, atomic force, and virial. As shown in Fig. S1 (see the supplementary material), the RMSEs of energy, force, and virial converge well with increasing training step.

We used the same hyperparameters for these two NEP models. After extensive testing, the selected hyperparameters were determined as follows: the radial and angular cutoffs are $r c R=8$ Å and $r c A=4Å$⁠, respectively. The number of radial and angular descriptor components is $n max R+1=9$ and $n max A+1=9$⁠, respectively. For angular parameters, we have $l max 3 b=4$ for three-body and $l max 4 b=2$ for four-body terms, respectively. The number of neurons in the hidden layer of the neural network is 50. The size of the population is $N pop=50$⁠, and the number of generations is $N gen=5× 10 5$⁠. The weights of energy, force, and virial RMSEs in the loss function were set to 1.0, 1.0, and 0.1, respectively.

All MD simulations were conducted at the target temperature using the Nose–Hoover chain method⁴⁹ with a time step of 1.0 fs. Initially, the simulation was run for 0.1 fs in the isothermal–isobaric ensemble and subsequently in the isothermal ensemble for 1.0 ns to achieve equilibrium. Afterward, HNEMD simulations were performed in the isothermal ensemble for 10.0 ns to calculate the running LTCs. We employed $4×16×8$ and $9×5×5$ supercells for $β− Ga 2 O 3$ and $κ− Ga 2 O 3$⁠, respectively, containing 10 240 and 9000 atoms. The simulation sizes were validated to be sufficiently large to eliminate finite-size effects. For each case, the predicted LTC was calculated as the average of five independent simulations, and the corresponding standard error was also estimated.

Figures 2(a)–2(c) show a comparison between the predicted energy, force, and virial values obtained from the NEP models and the corresponding DFT reference values for both the training and test sets of the two phases. In all scenarios, the RMSEs for energy, force, and virial are below 0.5 meV/atom, 40 meV/Å, and 5 meV/atom, respectively. These results demonstrate a very high accuracy of our NEP models. Additionally, in Fig. S2 of the supplementary material, we present the energy changes in $β− Ga 2 O 3$ and $κ− Ga 2 O 3$ by employing single atom perturbation, as predicted by NEP and DFT. We randomly select an oxygen atom, move it along the [100] direction, and calculate the change in total energy. The predictions made by NEP models closely align with the results from DFT calculations, further demonstrating the precision of our developed NEP models. In addition to their high accuracy, our NEPs are remarkably efficient. For a system comprising 10 000 atoms, our NEPs can attain a computational speed of approximately $1.5× 10 06$ atom-step per second in the GPUMD package using a single GeForce RTX 3090 GPU card. Compared with the BTE-ALD approach using ShengBTE⁵⁰ package with force constants obtained from DFT calculations,⁵¹ our GPUMD + NEP approach is much more computationally efficient (see Fig. S2 of the supplementary material), which enables us to perform extensive simulations to study the LTCs of the $Ga 2 O 3$ crystals.

To assess the reliability of our NEP models in capturing the phonon transport properties of $β− Ga 2 O 3$ and $κ− Ga 2 O 3$⁠, we compare the calculated phonon dispersions using both NEP and DFT methods, as shown in Fig. 3. It can be seen that for both phases, the acoustic branches predicted by NEP and DFT are very close, while the optical branches show slight deviations, especially at high frequencies. It should be noted that high-frequency optical phonons (>10 THz) contribute minimally to LTC,¹⁸ which is also confirmed by subsequent calculations of the spectrally decomposed LTC (see Fig. 6). Thus, we believe that the NEP models can reliably predict the LTC of the two phases of $Ga 2 O 3$⁠.

After confirming the accuracy of the NEP models, we apply them to calculate LTCs using the HNEMD method. Figure 4 presents the LTC for both $β$- and $κ− Ga 2 O 3$ using the HNEMD method at 300 K. Notably, Fig. 1(a) reveals that $β− Ga 2 O 3$ possesses a monoclinic crystal structure, suggesting the potential for non-zero values in the off-diagonal elements in the LTC tensor. Fortunately, the HNEMD approach [see Eq. (1)] allows us to fully determine the LTC tensor as listed in Table I, by applying the external driving force along three orthogonal axes sequentially. As anticipated, the $β$ phase exhibits minor couplings between the $x$ and $z$ directions. In contrast, the off-diagonal elements of the LTC tensor for the $κ$ phase are zero, attributed to its orthorhombic crystal structure.

For the $β$ phase, the LTC is highest along the [010] direction with $19.9±0.2$ W/(m K), followed by [001] with $12.6±0.2$ W/(m K) [see the blue dotted line in Fig. 5(e)] and [100] at $10.3±0.2$ W/(m K). As shown in Table II, our HNEMD results are consistent with previous experimental results^52,53 using time-domain thermoreflectance and theoretical predictions based on the BTE-ALD method^{18–20,51,54} or the EMD method.^15,27 This further demonstrates the reliability of our HNEMD approach based on machine-learned NEP in characterizing the LTC of $Ga 2 O 3$ crystals.

Furthermore, we predict that the LTC for $κ− Ga 2 O 3$ to be $4.5±0.1$⁠, $3.9±0.1$⁠, and $4.0±0.1$ W/(m K) along the [100], [010], and [001] directions, respectively. The results are consistent with those of a previous work which combines BTE with ALD to study the thermodynamic properties of $Ga 2 O 3$⁠.⁵¹ These values are about one-fifth to half of those of $β− Ga 2 O 3$⁠. While the $β$ phase displays a significant LTC anisotropy with an anisotropy index of 1.94 [the maximum-to-minimum ratio of the LTC, see spatial distribution in Fig. 5(b)], the $κ$ phase exhibits a nearly isotropic LTC pattern with an anisotropy index of 1.15 [see Fig. 5(c)].

Figure 6 presents the spectrally decomposed LTC of $β$- $Ga 2 O 3$ and $κ− Ga 2 O 3$ calculated from Eq. (5). Taken as a whole, only phonon modes with frequencies in the 0–10 THz range are really involved in thermal transport. It also can be found that, although in both phases the low-frequency acoustic phonons contribute significantly to the LTC, the $β$ phase has a much broader distribution of phonon frequencies contributing to its LTC than the $κ$ phase. This is also related to the phonon dispersions as shown in Fig. 3. The branches of the $κ$ phase with a phonon frequency in the range of 5–10 THz are much flatter than those of the $β$ phase, leading to much lower phonon group velocities. Moreover, in the $κ$ phase, we observe significant overlaps between multiple bands for phonon frequencies greater than 3 THz, indicating the presence of multiple scattering channels, a phenomenon previously observed in Violet phosphorene.²⁶ 

Furthermore, Fig. 7 shows the temperature-dependent LTC for $β$- and $κ− Ga 2 O 3$⁠. In the temperature range of 200–600 K, the LTCs of both phases decrease with increasing temperature. The LTC of the $β$ phase exhibits a temperature dependence slightly stronger than $∼ T − 1$⁠, as typical for systems dominated by three-phonon scattering processes.⁵⁵ However, the LTC of $κ− Ga 2 O 3$ shows a clearly weaker temperature dependence, which is $∼ T − 0.5$ along the [100] direction and $∼ T − 0.7$ along the [010] and [001] directions, suggesting a significant impact of higher-order anharmonic phonon scatterings attributable to its complex crystal structures.⁵⁶ 

In summary, we have developed machine-learned NEP models trained against quantum-mechanical DFT data for the $β$ and $κ$ phases of $Ga 2 O 3$⁠, which have been demonstrated to be accurate and efficient in predicting energy, atomic forces, virial, and phonon dispersions in both phases. Based on large-scale HNEMD simulations, we reached a consistent prediction with previous experimental measurements for the $β$ phase and predicted the LTC of the $κ$ phase. We found that the $κ$ phase has a much lower LTC than the $β$ phase, due to its phonon frequency contributions being limited to below 5 THz compared to the 0–10 THz range of the $β$ phase. Furthermore, based on the tensor analysis, the $κ$ phase exhibits an almost isotropic spatial distribution with an anisotropy index of 1.15 at 300 K, in contrast to the pronounced anisotropy index of 1.94 for the $β$ phase. We also examined the temperature dependence of the LTC in the two phases and found that the LTC of the $β$ phase follows a temperature dependence slightly stronger than $∼ T − 1$⁠, whereas the $κ$ phase shows a weaker temperature dependence from $∼ T − 0.5$ to $∼ T − 0.7$⁠, indicating a significant effect of high-order anharmonicity as in low-LTC materials. Our work demonstrates that the machine-learned NEP-driven HNEMD simulations can reliably and effectively characterize phonon thermal transport properties for complex crystals such as $κ− Ga 2 O 3$⁠, and thus we expect that this approach can be used to explore the LTCs of other phases of $Ga 2 O 3$ as well.

See the supplementary material for the evolution of various loss functions during the training process, a comparison of computational efficiency between the GPUMD + NEP potential and the DFT + BTE-ALD method, and the validation of energy changes with single atom perturbation against DFT calculations.

This work was supported by the Key-Area Research and Development Program of Guangdong Province (Grant No. 2020B010169002), the Guangdong Special Support Program (Grant No. 2021TQ06C953), and the Science and Technology Planning Project of Shenzhen Municipality (Grant No. JCYJ20190806142614541).

The authors have no conflicts to disclose.

Xiaonan Wang: Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (lead). Jinfeng Yang: Investigation (supporting); Methodology (supporting); Software (supporting); Writing – review & editing (supporting). Penghua Ying: Investigation (supporting); Methodology (equal); Software (equal); Validation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Zheyong Fan: Methodology (supporting); Software (supporting); Writing – review & editing (equal). Jin Zhang: Writing – review & editing (supporting). Huarui Sun: Conceptualization (lead); Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal).

Complete input and output files for the NEP training of $β− Ga 2 O 3$ and $κ− Ga 2 O 3$ crystals are freely available at https://gitlab.com/brucefan1983/nep-data. The source code and documentation for GPUMD are available at https://github.com/brucefan1983/GPUMD and https://gpumd.org, respectively.