Machine-learned atomic cluster expansion potentials for fast and quantum-accurate thermal simulations of wurtzite AlN Open Access

Wurtzite aluminum nitride (w-AlN) emerges as a promising semiconductor, distinguished by various exceptional characteristics. These include an ultrawide bandgap^1–3 (∼6.1 eV), a large critical electric field^1,3 (∼15 MV/cm), a high sound velocity⁴ (∼11 km/s), large piezoelectric coefficients,⁵ a relatively high thermal conductivity^6,7 (⁠ $κ∼300 W / m K$⁠), and a lattice similar to other semiconductors such as w-GaN. The large critical electric field stemming from the ultrawide bandgap results in Baliga's figure of merit (FOM) and Johnson's FOM of w-AlN significantly surpassing those of w-GaN or β-Ga₂O₃,^1,3 thereby establishing w-AlN as a remarkable candidate for novel high-power or radio frequency (RF) electronics.^8,9 Meanwhile, the ultrawide bandgap of w-AlN facilitates developments in deep-ultraviolet photonics.^3,8,10,11 The high sound velocity and piezoelectric performance render w-AlN suitable for fabricating microelectromechanical system (MEMS)-based resonators and filters, which are extensively applied in 5G communications.^12–15 In addition, the thermal conductivity of w-AlN is considered satisfactory; hence, w-AlN sometimes serves as high-κ substrates for high-power devices to improve the heat dissipation performance.^9,16,17 As depicted in Fig. 1, the bandgap and thermal conductivity of w-AlN are compared to other representative materials,² further exhibiting the significant promise of w-AlN.

Though the inherent characteristics of w-AlN are excellent, further studies on its properties remain necessary for developing next-generation electronics, especially in the aspects of (i) clarifying how the complex environment affects the physical properties, (ii) tuning the physical properties on demand, and (iii) analyzing the growth processes. One specific issue of interest pertains to the effects of lattice strains on the thermal properties of w-AlN. Since w-AlN is extensively used as the nucleation layer within the GaN high electron mobility transistors (HEMTs) to buffer lattice mismatch between w-GaN epilayers and substrates,^22–25 the residual stress within AlN is inevitable. However, the correlations between lattice strains with the thermal conductivity and phonon bands of w-AlN remain vague, which may affect the heat dissipation and reliability of corresponding devices.^26,27

In addition to experimental approaches, atomistic simulations act as another avenue for gaining insights into the physical properties of novel materials,²⁸ which is traditionally represented by two techniques,²⁹ i.e., the first-principle calculations based on density functional theory (DFT) and the molecular dynamics (MD) simulations based on empirical potentials. Nevertheless, high computational cost limits the DFT methods for modeling transport properties, while the MD simulations based on simple empirical potentials are less accurate than DFT.³⁰ For w-AlN, several empirical potentials have been proposed, including the Stillinger–Weber (S–W),³¹ Tersoff,³² Vashishta,³³ and COMB3³⁴ models. However, each empirical potential generates divergent lattice parameters or phonon dispersions from those of DFT.^33–35 When predicting thermal conductivity, it is fundamental to accurately describe both harmonic and anharmonic interactions of phonons.^36–38 This places a heightened demand on the accuracy of interatomic potentials for w-AlN.

In recent years, machine learning (ML) interatomic potentials have attracted significant attention by effectively balancing computational efficiency with accuracy. A wealth of literature has shown that a well-built ML potential trained with the DFT reference data can provide an unbiased representation of potential energy surfaces and simultaneously exhibit strong transferability.^36,39,40 More importantly, the linear behavior in the computational cost of ML potentials enables them with much higher efficiency and scalability than DFT methods. Until now, several ML potential models have been proposed, such as the neural network potential (NNP),^41,42 Gaussian approximation potential (GAP),^40,43,44 spectral neighbor analysis potential (SNAP),^45,46 deep potential (DP),^47,48 moment tensor potential (MTP),^49,50 atomic cluster expansion (ACE) potential,^51,52 neural equivariant interatomic potential (NequIP),⁵³ Allegro,⁵⁴ and MACE.^55,56 The ACE potential is one of the most computationally efficient and quantum-accurate models available⁵² and is also suitable for performing large-scale simulations on CPU platforms. Hence, ACE is chosen in this work.

Here, we introduce a machine-learned ACE potential^51,52,57 for w-AlN, aiming to facilitate the atomistic simulations of its thermal properties and gain insights into how tuning factors (such as lattice strains) influencing the phonon pictures. The remainder of this paper is organized as follows. In Sec. II, we concisely introduce the ACE methodology and the construction of the training database for w-AlN. In Sec. III, we comprehensively demonstrate the accuracy of our ACE potential in predicting various thermal and mechanical properties of w-AlN by comparison with either DFT calculations or experiments. Then, our ACE potential is applied to unravel the correlations between thermal conductivities and biaxial strains of w-AlN. Essential conclusions of this study are presented in Sec. IV.

The coefficients $c$ are free model parameters that can be optimized during fitting. Further information about the ACE architecture can be found in the literature.^51,52,57 In this work, we employ the software package Pacemaker⁵⁸ for the parametrization of the ACE potential. The final ACE model associated with the detailed hyperparameters is freely available in the supplementary material.

Structures in the training database are obtained from MD trajectories. Here, the empirical S-W potential³¹ is adopted to carry out the MD simulations. This approach allows for sampling over extended time scales (∼ns) to ensure structural diversity, while it bypasses the computationally intensive ab initio MD and makes the sampling more efficient. In detail, the initial structure consists of a 3 × 3 × 3 supercell containing 108 atoms. The cell is then expanded or compressed with a scaling factor ranging from 0.95 to 1.05 on each lattice constant. Next, a series of MD simulations for each cell are performed with the canonical (NVT) ensemble at temperatures of 100, 500, and 1000 K. For each temperature, a 1-ns trajectory is produced, from which structures are sampled at uncorrelated intervals of 150 ps. Finally, a total of 13 608 local atomic environments are collected from all 21 MD trajectories.

All generated structures are subsequently subjected to single-point DFT calculations to obtain well-converged reference energies and forces for training. The DFT calculations are performed with the Vienna ab initio simulation package (VASP).⁵⁹ Exchange and correlation are treated by using the PBEsol functional⁶⁰ with a projector-augmented wave method.⁶¹ Moreover, we adopt Gaussian smearing of 0.05 eV width to electronic levels, a 600-eV cutoff for plane wave expansions, and a maximum spacing of 0.2 Å⁻¹ for meshing the reciprocal space. The total energy is attained with a convergence criterion of less than 10⁻⁶ eV in the self-consistent electronic iterations.

The accuracy of our ACE potential is validated through a comparison with DFT-predicted energies and forces. The training and testing datasets comprise 36 612 and 4212 atomistic force components, respectively. In Fig. 2, we present a comparison of the total energies and atomic forces predicted by our potential with those derived from DFT. Notably, our ACE potential effectively reproduces the total energies, showcasing a remarkably low root-mean-squared error (RMSE) of 0.13 meV/atom for the testing datasets. Furthermore, the interatomic forces in the testing datasets are accurately predicted, with a relatively low RMSE of 5.01 meV/Å. The accuracy of both energy and force predictions achieved with our ACE potential is comparable to those of other reported works,^{40,43,52,62,63} which demonstrates that our ACE potential serves as a robust representative of the DFT potential energy surface.

Then, we assess the capability of our ACE potential to model the physical properties of w-AlN. Our evaluation begins by employing the ACE to predict lattice parameters, a fundamental property that significantly influences various intrinsic characteristics of a material. The lattice parameters produced by our ACE potential are as follows: $a=b=3.11558Å$⁠, $c=4.98513Å$⁠, and $γ= 120 °$⁠. These values are in excellent agreement with both our DFT calculations and experimental results⁶⁴ (Table I). For comparison, lattice parameters are also predicted using the S-W potential³¹ here as well as using the Tersoff potential elsewhere.³⁵ However, it is noteworthy that the relative errors in these predictions are considerably larger than those obtained with ACE.

The ACE potential is subsequently employed to predict the volumetric specific heat capacity (⁠ $C V$⁠), coefficients of thermal expansion, and the bulk modulus of w-AlN. $C V$ is calculated using the Phonopy package^65–67. Figure 3(a) summarizes the $C V$ obtained from our ACE potential, the DFT calculations, and the experimental measurements.⁶⁸ Obviously, the data generated by our ACE potential exhibit strong agreement with both DFT and experiments.

Moreover, a mechanical property, namely, the bulk modulus (⁠ $B T$⁠) of w-AlN, is also determined by the ACE potential. As shown in Fig. 3(c), the ACE potential quantitatively reproduces the $B T$ as DFT calculations. Both ACE and DFT predict a subtle decrease in the modulus as the temperature increases from 0 to 1000 K. In summary, it is obvious that our ACE potential is capable of accurately describing both thermal and mechanical properties.

The prediction of the phonon dispersions of a material is another crucial metric for the quality of a potential to describe lattice dynamics. We first calculate the second-order harmonic and third-order anharmonic interatomic force constants (IFCs) through the finite displacement method.⁶⁹ By combining the Phonopy package with the second-order IFCs calculated from our ACE and DFT, the phonon dispersions of w-AlN at 0 K are determined, as illustrated in Fig. 3(d). As a characteristic of the ionic crystal, the splitting of LO–TO phonons at the Γ-point is observed, which is attributed to the long-range Coulomb interactions.³⁶ Since our ACE model does not encompass the Coulomb interactions, we have incorporated the non-analytical correction²⁶ into the dynamical matrix to resolve the splitting of LO–TO phonons at the Γ point. Our results show that the ACE model accurately predicts the phonon frequencies at almost of all high-symmetry points and accurately captures the dispersion behavior of each phonon branch. Meanwhile, the phonon density of states (DOSs) calculated by ACE is almost identical to the DFT results.

In fact, although the absence of Coulomb terms introduces some errors in force calculations, it will be shown later that the accuracy of our potential remains sufficiently high to yield thermal conductivity predictions. If the long-range Coulomb interactions become crucial for other calculations, we will properly incorporate the Coulomb interactions with fixed partial charges before fitting the ACE model. This adjustment will automatically resolve the LO–TO splitting issue and improve the accuracy of force calculations.³⁶ Similar treatments for GaN can be found in the literature.⁴³ 

Herein, we utilize the second- and third-order IFCs to predict the thermal conductivity (⁠ $κ$⁠) of w-AlN. The Wigner transport equation (WTE)^70–72 is solved by the direct solution in Phono3py^66,73,74 package, with a 25 × 25 × 25 mesh for sampling the first Brillouin zone over temperatures within 100–600 K. Due to its wurtzite lattice, it is expected that the thermal conductivities of w-AlN exhibit approximate isotropy along both in-plane and cross-plane directions.²⁷ 

We also conduct ShengBTE⁷⁵ calculations including the four-phonon processes (4ph)⁷⁶ for comparison with the WTE. Since the computational cost is extremely huge with 4ph involved, we resort to the sampling-accelerated method⁷⁷ with the settings “num_sample_process_4ph = 1E5” and “num_sample_process_4ph_phase_space = 1E5.” Meanwhile, a 16 × 16 × 16 mesh of the first Brillouin zone is set, and the value of “scalebroad” is taken as 0.1. The cutoff distances of the third- and fourth-order IFCs are selected to be the seventh nearest and the second nearest atom, respectively. These mentioned parameters are large enough to make it converge to the rigorous calculations,^75–78 though resulting in a negligible uncertainty due to random sampling.

As illustrated in Fig. 4(a), the ACE potential calculates the thermal conductivities of w-AlN along both in-plane and cross-plane directions. Though the physical pictures behind WTE and 4ph-scattering are different, their predictions are close indeed. The thermal conductivity predicted by ACE shows overall good agreements with our DFT calculations and experimental measurements (based on the three-sensor 2ω method proposed in our recent works).^7,22,80 Besides, the ACE-predicted values approximate the literature experiments^2,79 above room temperature as well, though the ACE overestimates thermal conductivities beneath 250 K owing to neglecting the phonon-defect scattering. Several theoretical and experimental studies^27,79,81,82 have discussed this issue, and researchers generally attribute the discrepancy between theoretical predictions and cryogenic experiments to the impacts of phonon-defect (e.g., point defects, dislocations, grain boundaries) scattering. Since the phonon-defect scattering rates are independent of temperature in principle,²⁷ they will dominate the phonon transport at low temperatures where the anharmonicity (normal and Umklapp processes) is weak.

The literature experiments depicted here are all based on an isotropic assumption, whereas our three-sensor 2ω method enables the direct derivation of thermal conductivities along different directions.⁷ The brief introduction on the three-sensor 2ω method can be found in the supplementary material.

Furthermore, we conduct an equilibrium molecular dynamics (EMD) simulation based on the ACE potential to predict the w-AlN thermal conductivity. The large quantum effects at low temperatures significantly influence the MD-calculated thermal conductivity,^{28,81,83–86} since the Debye temperature of w-AlN is quite high (⁠ $Θ D∼1000 K$⁠).⁸⁷ Hence, we only calculate the thermal conductivity at 1100 K, where the quantum effects could be reasonably neglected.⁸⁴ A bulk w-AlN system containing 4000 atoms is set as the initial configuration for EMD simulations conducted in the LAMMPS.⁸⁸ The length of each crystallographic direction reaching 10 unit cells achieves convergence for calculating the thermal conductivity of w-AlN.^89,90

Periodic boundary conditions are applied to all three directions to mimic the infinite size of structures. A total of 4.5 ns EMD simulations based on the trained ACE potential are performed with a time step of 1 fs. After equilibrating the system in the isothermal–isobaric (NPT) ensemble for 500 ps, the system is switched to the microcanonical (NVE) ensemble for 4 ns to collect heat flux data. To mitigate the impacts of statistical uncertainties and errors, we conduct ten independent EMD simulations with different initial velocity distributions to obtain the averaged thermal conductivities at 1100 K. The thermal conductivity here is considered isotropic along the in-plane directions, so we average the thermal conductivity values along the x, y directions as the final in-plane thermal conductivity. The EMD results are $κ in MD=51.0±5.8 W / m K$ and $κ cr MD=48.4±6.3 W / m K$⁠, which are comparable to those of BTE calculations, as discussed in the supplementary material. On the AMD EPYC^TM 7452 CPU platform, an average calculation efficiency of the ACE-based EMD simulations reaches 0.26 ms/MD-step/atom when calculating on a single thread (using GCC 9.1.0 compiler and LAMMPS 23-Jun-2022-Update3), which manifests the capability of the ACE potential for large-scale molecular dynamics.⁵² 

More detailed phonon transport characteristics of w-AlN are also calculated from the ACE potential by Phono3py. The accumulated thermal conductivity as a function of phonon frequency at 300 K is presented in Fig. 4(b), which further verifies the accuracy of the ACE potential compared to the DFT results. Obviously, the thermal conductivity is primarily contributed by acoustic phonon branches with frequencies below the “phonon bandgap,” whose contributions exceed 98.5% for both in- and cross-plane values. In addition, we also calculate the accumulated thermal conductivity as a function of phonon mean free path (MFP) based on the single-mode relaxation time approximation, as detailed in the supplementary material. The discussions on MFP imply that the size effect of w-AlN film's thermal conductivity is crucial for practical applications, which limits the heat dissipation performance of the corresponding electronic devices.^{22,23,25,91,92}

As discussed in Sec. I, there is an inevitable residual strain (stress) inside w-AlN in its practical applications, especially when w-AlN serves as the transition layer for a GaN HEMT by forming the GaN/AlN/Substrate heterojunction. Owing to the mismatches in lattice and thermal expansion between three materials, residual strain and lattice defects inside the w-AlN are general significant. Therefore, based on the trained ACE potential of w-AlN, we proceed to study the strain effects on lattice thermal conductivity of w-AlN.

From a classical and intuitional perspective, compressive strain is expected to increase the thermal conductivity by augmenting the elasticity modulus and acoustic velocity.²⁶ The results, shown in Fig. 5, reveal that both in- and cross-plane thermal conductivity of w-AlN decrease remarkably under the +4% biaxial strain state (tensile) at room temperature, with the average thermal conductivity decreasing by 40%. Conversely, under the −4% biaxial strain state (compressive), the average thermal conductivity increases by 30%. In the temperature range of 200–400 K, lattice thermal conductivity decreases significantly; however, the decrease/increase trend depending on biaxial strain is nearly constant at different temperatures [Fig. 5(b)], implying that the influences of lattice strain and temperature on thermal conductivity should be independent. Meanwhile, the temperature dependence of thermal conductivity follows a similar trend under different biaxial strain states.

Figures 6(a) and 6(b) illustrate the phonon harmonic properties under different strain states, while Figs. 6(c) and 6(d) depict the anharmonic features. Volumetric specific heat of phonon modes reflects the energy level of a crystal system, varying with the phonon dispersions and volume of the unit cell correspondingly. As shown in Fig. 6(a), slight decreases of mode specific heat with tensile biaxial strains occur, implying a positive correlation to the change of thermal conductivity. In Fig. 6(b), phonon group velocity gradually increases from the tensile to compressive strain states, which is consistent with the increased thermal conductivity. Figure 6(c) shows that phonon relaxation time decreases under the tensile strain state and increases under the compressive strain state, which is consistent with the thermal conductivity variations as well. According to the results of phonon DOSs in the supplementary material, the phonon bandgaps are 0.65, 1.35, and 2.40 THz under tensile, free, and compressive strains, respectively. This determines the variations of relaxation time in principle, since a smaller (larger) phonon bandgap will enable more (less) available three-phonon scattering channels, thereby enhancing (suppressing) the three-phonon anharmonic scattering processes.^26,95

The Grüneisen parameter, an indicator of lattice anharmonicity,^26,27 is further analyzed here to back up the influences of strains on anharmonic scattering processes [Fig. 6(d)]. Obvious increases (decreases) in the Grüneisen parameters of acoustic phonon branches are observed under the tensile (compressive) strain state, in line with the variations of phonon bandgaps. More discussions on how and why lattice strains affect the phonon anharmonicity can be found in the supplementary material. Consequently, the consistent variations in mode specific heat, group velocity, relaxation time, and anharmonicity lead to significant changes in the lattice thermal conductivity²⁶ of w-AlN under different strain states.

The highlighted impacts of lattice strains on phonon properties should be favorable for tuning the heat dissipation performance of corresponding w-AlN-based electronic devices. By carefully engineering the residual strains within heteroepitaxial structures via annealing^96,97 or selecting substrates with specific lattice structures,^26,98 it is viable to optimize the lattice thermal properties on demand, facilitating near-junction thermal optimization of semiconductor devices.^91,99

We have developed an ACE potential based on machine learning for atomistic simulations of monocrystalline w-AlN. Our ACE potential exhibits remarkable accuracy in reproducing the DFT potential energy surface of w-AlN, achieving an energy RMSE of ∼0.13 meV/atom and a force RMSE of ∼5.01 meV/Å for both Al and N atoms. Subsequently, the predictive power of ACE is demonstrated across a variety of properties of w-AlN, including ground-state lattice parameters, specific heat capacity, coefficients of thermal expansion, bulk modulus, phonon dispersions, and thermal conductivity. All these results show excellent agreement with the DFT calculations and experimental results, demonstrating that the ACE model sufficiently describes both harmonic and anharmonic phonon properties.

The lattice strain is proven as a significant tuning factor for thermal design of heteroepitaxial electronic devices. As a practical application of the ACE potential, we perform lattice dynamics simulations to unravel the effects of biaxial strains on thermal conductivity of w-AlN. The results indicate that a 4% biaxial tensile (compressive) strain approximately causes a 40% decrease (30% increase) in the thermal conductivity of w-AlN, while the influences of lattice strain and temperature on thermal conductivity appear to be independent. The investigations into phonon pictures under different strains reveal that the consistent variations in phonon mode heat capacity, group velocities, and relaxation times predominantly contribute to the variation of thermal conductivity, while all these factors stem from the variations in phonon band structures. Thus, it is feasible to facilitate near-junction thermal optimization of devices via strain engineering on phonon bands. The findings here are beneficial for the development of next-generation electronic devices.

See the supplementary material for brief introductions to the experiments, additional discussions on phonon properties, and the ACE potential files developed for simulations in both LAMMPS and ASE calculators.

This work was financially supported by the National Natural Science Foundation of China (NNSFC) (Grant Nos. U20A20301, 52327809, 52250273, and 51825601).

The authors have no conflicts to disclose.

Guang Yang and Yuan-Bin Liu contributed equally to this work.

Guang Yang: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Software (equal); Validation (lead); Visualization (lead); Writing – original draft (lead). Yuan-Bin Liu: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (supporting); Methodology (equal); Resources (supporting); Software (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Lei Yang: Methodology (supporting); Software (supporting). Bing-Yang Cao: Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The ACE potential is attached in the supplementary material, which can be used with the ML-PACE package on https://github.com/ICAMS/lammps-user-pace. The data that support the findings of this study are available from the corresponding author upon reasonable request.