Large-Scale MD Simulation of GaN/AlN Heteroepitaxial Thin-Film Cracking Using LightPFP

Introduction

Fig. 1. Motivation: thermal-cooling-induced cracking in a GaN/AlN/substrate heterostructure.

GaN-on-AlN thin film systems are important for next generation power electronics, communication devices, and UV LEDs. Cracks occurring after epitaxial growth is an eminent challenge in device fabrication (Fig. 1).

The fundamental questions are: where does cracking initiate, and does the AlN substrate orientation matter? This problem could not be addressed with confidence by existing simulation approaches. Classical Stillinger-Weber (SW) potentials are fast but their intermolecular interaction especially at the interface are inaccurate. PFP is a major improvement in terms of accuracy but is limited to tens of thousands of atoms.

Cracking at the interface requires a larger system to not only simulate crack nucleation but also crack propagation. We trained a LightPFP (LPFP) model for a GaN/AlN system using PFP generated data coupled with active learning.

LPFP can run simulations up to hundreds of thousands of atoms. This work demonstrates that LPFP quantitatively reproduces PFP-quality crack initiation and propagation behavior at the length scales required for engineering-relevant simulations.

We also investigated two AlN substrate orientations and found that they produce fundamentally different cracking mechanisms.

The orientation dependence has direct implications for crack initiation in polycrystalline GaN films used in commercial devices.

Computational details

GaN on AlN interface was initially equilibrated for two AlN orientations: polar c-plane (0001) and non-polar a-plane (11-20) (Fig. 2)¹. System sizes of 1x1x1 (supercell), 2x2x1, 3x3x1, and 5x5x1 were evaluated (Fig. 3). Classical Stillinger-Weber potential (SW)², PFP, and LightPFP (LPFP) were compared.

The LPFP model was trained using PFP (v8.0.0, r2scan) as the data source. Initial dataset: 300 structures from bulk wurtzite GaN, AlN, and a/c-plane interface structures. Active learning iterations collected around 200 additional structures. Periodic boundaries were in the x, y direction and z was a free surface.

Cracking simulations used LAMMPS³ NPT with fix deform (x erate 0.002-0.2 ps⁻¹) at 300 K for system sizes 7,896-197,400 atoms.

Polyhedral Template Matching in OVITO⁴ was used to calculate the nonhexagonal bond ratio (relative disorder fraction) to quantitatively analyze the onset of cracking.

Figure 2. Atomic structures of the two AlN orientations studied: a-plane (left) and c-plane (right).

Figure 3. GaN/AlN supercells (7,896 atoms to 197,400 atoms) for size-dependence analysis of the system.

Results & Discussion

The crack initiation of LPFP shows good agreement with PFP unlike SW (Fig. 4). The fluctuation in the low strain elastic regime for SW can be attributed to low stability of wurtzite ordering. Cracking starts at a similar strain for both force fields at all system sizes, which signifies size convergence and rules out out finite-size artifacts for LPFP. Though crack propagation characteristics are different between LPFP and PFP, both eventually plateau at the same value.

The system size difference for cracking was smaller for LPFP than SW, which maybe due to the initial disorder fluctuation (Fig. 5). Within a specified region, the relative disorder fraction is the number of non-hexagonal atoms at each strain normalized by the initial number of non-hexagonal atoms.

Fig. 4. Relative Disorder Fraction vs. strain for LPFP, PFP, and SW (a-plane, 1×1×1); inset shows SW's spurious elastic-regime fluctuations.

Fig. 5. System-size effect on Relative Disorder Fraction: LPFP (left) scales consistently, while SW (right) shows erratic, size-dependent behavior.

The a-plane system can sustain more strain than the c-plane, which is essential information in preparing the GaN epitaxial layer (Fig. 6). This difference originates at the GaN/AlN interface for the c-plane system. Though the a-plane has numerous point defects at the onset in the GaN layer there are no line dislocations in the interface and therefore shows higher strain tolerance (Fig. 7).

The line defects in the c-plane act as growth sites, which thicken into planes gradually with strain, and cracks clearly propagate vertically. On the other hand, a threshold strain is necessary for crack nucleation to occur for the a-plane system and the source is likely from the AlN layer rather than the interface. The plateaus after the crack occurring also signifies proper crack propagation characteristics unlike a-plane SW. Note the fluctuation of the disorder fraction in the elastic regime of c-plane SW. Though the c-plane initial interfacial dislocation is identical for both methods, SW has low crystallinity with increasing strain, which may be from imprecise energy landscape mapping. This is clear evidence of LPFP being more accurate than SW.

Fig. 6. Orientation dependence of fracture onset: LPFP vs. SW for c-plane and a-plane (300 K, 0.02 ps⁻¹).

Fig. 7. Crack nucleation and propagation under LPFP for c-plane (left) and a-plane (right); fully hexagonal atoms omitted.

Calculation conditions

Items	Details
Number of atoms	7,896 (1x1x1) to 197,400 (5x5x1) atoms; c-plane 189,039; a-plane 197,400
Elements	Ga, N, Al (GaN on AlN, wurtzite structure)
Potential	LightPFP, PFP (v8.0.0 r2scan), SW
MD	LAMMPS NPT: 0Pa 300K p p f (periodic XY, free Z), reflective z-walls deform x strain rate: 0.2, 0.02, 0.002 ps^-1

References

[1] R. Li, G. Wu, K. Liang, S. Wang, X. Sun, X. Han, L. Xue, H. Li, S. Liu, Effect of AlN substrate orientation on GaN film quality, Comp. Mat. Sci. 202, 110991 (2022).

[2] H. Li, H. Xiang, H. Huang, Z. Zeng, X. Peng, Interface structure deformation mechanisms of AlN/GaN multilayers. Cer. Int. 46, 11556 (2020)．

[3] A.P. Thompson et al., LAMMPS — A flexible simulation tool for particle-based materials modeling. Comp. Phys. Comm. 271, 108171 (2022).

[4] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO – the Open Visualization Tool, Modelling Simul. Mater. Sci. Eng. 18, 015012 (2010).