Characterizing Defect Dynamics in Silicon Carbide Using Symmetry-Adapted Collective Variables and Machine Learning Interatomic Potentials

TL;DR

This paper tackles the challenge of simulating defect dynamics in SiC divacancies, which are hindered by large activation barriers, by developing an E(3)-equivariant Allegro MLIP trained through a six-round active learning workflow that couples symmetry-adapted PINES collective variables with enhanced sampling and a stabilizing prior. The resulting MLIP achieves ab initio accuracy and high efficiency (≈6 ns/day) in 216-atom SiC systems, accurately reproduces defect transition free energy barriers within thermal fluctuations (≈$k_B T \,\approx\,0.13$ eV at 1500 K), and transfers to larger multi-defect configurations. Temperature-dependent analyses using Markov State Models reveal five defect macrostates, with the divacancy showing maximum thermodynamic stability near ~1625 K, providing atomistic insight into annealing protocols for defect stabilization. The study demonstrates a generalizable framework for data-efficient MLIP development in condensed matter, offering a powerful tool for defect engineering in SiC and highlighting avenues for extending to charge/spin-aware potentials and cross-polytype transferability.

Abstract

Silicon carbide (SiC) divacancies are attractive candidates for spin defect qubits possessing long coherence times and optical addressability. The high activation barriers associated with SiC defect formation and motion pose challenges for their study by first-principles molecular dynamics. In this work, we develop and deploy machine learning interatomic potentials (MLIPs) to accelerate defect dynamics simulations while retaining ab initio accuracy. We employ an active learning strategy comprising symmetry-adapted collective variable discovery and enhanced sampling to compile configurationally diverse training data, calculation of energies and forces using density functional theory (DFT), and training of an E(3)-equivariant MLIP based on the Allegro model. The trained MLIP reproduces DFT-level accuracy in defect transition activation free energy barriers, enables the efficient and stable simulation of multi-defect 216-atom supercells, and permits an analysis of the temperature dependence of defect thermodynamic stability and formation/annihilation kinetics to propose an optimal annealing temperature to maximally stabilize VV divacancies.

Figures (8)

• Figure 1: Vacancy defects in silicon carbide (SiC) crystals. (A) Schematic illustration of the divacancy formation process in SiC. Vacancies are typically created by ion implantation or electron irradiation followed by temperature annealing to promote high defect mobility and the formation of stable divacancies. (B) A number of vacancy defects commonly appear and interconvert during the annealing process including divacancies (VV), antisite vacancies (V$_{C}$V$_{Si}$), carbon vacancies (V$_{C}$), and silicon vacancies (V$_{Si}$).
• Figure 2: Active learning protocol for training an Allegro machine learning interatomic potentials (MLIP) for SiC. The cyclic process includes (1) enhanced sampling using an autoencoder-based collective variable (CV), (2) configuration selection, (3) density functional theory (DFT) evaluation of selected configurations, (4) convergence assessment, (5) training and validation data appending and (6) force field retraining with the updated dataset.
• Figure 3: Accuracy of the trained Allegro MLIP after each round of active learning. Box plots of the mean average error (MAE) in the force predictions relative to the DFT ground truth over (A) the 300 training configurations and (B) 100 validation configurations in each round of the active learning cycle. The green horizontal line within the box represents the mean of the distribution, and the orange horizontal line represents the median. The limits of the box correspond to the first (Q1) and third (Q3) quartiles, and the whiskers extend to the minimum and maximum values within 1.5 times the interquartile range (IQR) from the quartiles. Outliers in the distribution are explicitly shown. The red line connecting the means indicates the trend. The dotted horizontal line represents the desired force MAE threshold of 0.1 eV/Å used to judge convergence of the model. (C,D) Analogous plots to panels (A,B), but illustrating the MAE in the per atom energy relative to the DFT ground truth. The dotted horizontal line represents an MAE in the energy of 0.005 eV/atom. Energy distribution of each round of (E) training and (F) validation configurations showing the diversity of configurations explored over each of the active learning rounds. The exploration of high-energy training configurations in Round 2 coincides with the spike in the MAE.
• Figure 4: Free energy landscapes for five defect transition processes estimated under the trained Allegro MLIP using ABF: (A) reorientation of a VV divacancy via motion of a Si atom, (B) migration of a V$_{C}$ monovacancy, (C) migration of a V$_{Si}$ monovacancy, (D) VV divacancy formation from a pair of monovacancies via C-atom hopping, and (E) VV divacancy formation from a pair of monovacancies via Si-atom hopping. In each case, the free energy landscape is projected into one dimensional projection CV as described by Zhang et al.Zhang2024.
• Figure 5: Force accuracy of trained Allegro MLIP in extrapolative predictions to (A) $V_C V_{Si} V_C$ and (B) $V_{Si}$V$_C$V$_{Si}$ trivacancy defect simulations and (C) multivacacancy defect simulations. The MAE of the MLIP force predictions relative to the DFT calculations are illustrated over the simulation time series. The black lines and shaded regions represent the MAE and the standard deviation of the per-atom force error, respectively, computed at 100 ps time intervals.