Virtual Expo 2026

Molecular Dynamics Simulation of Radiation Damage and Self Healing Property of Gallium Nitride (GaN)

Envision Piston

google meet link : https://meet.google.com/cjz-ycty-ace

1. Introduction

Gallium nitride is a wide bandgap semiconductor (E_g approximately 3.4 eV) that has attracted significant research attention for use in power electronics, high electron mobility transistors (HEMTs), light emitting diodes, and radiation detectors. In all of these applications, GaN can be exposed to energetic particle radiation, whether neutrons in a nuclear reactor, protons and heavy ions in space, or electrons in a particle detector. These energetic particles interact with the GaN lattice through nuclear collisions, transferring enough energy to displace atoms from their lattice sites and creating point defects.

The fundamental event in radiation damage is the primary knock on atom (PKA) cascade. A PKA is an atom that receives a large kinetic energy from an incoming particle and proceeds to collide with neighboring atoms, creating a branching chain of secondary and tertiary collisions called a displacement cascade. The cascade produces a mix of vacancies (empty lattice sites) and interstitials (atoms sitting between lattice sites), collectively known as Frenkel pairs.

What makes GaN particularly important to study is its demonstrated radiation tolerance. Unlike silicon or gallium arsenide, GaN shows substantial self healing behavior: many of the Frenkel pairs created during a cascade recombine spontaneously at room temperature through thermally activated migration. This recovery reduces the permanent defect concentration below what simpler models would predict. Understanding this process at the atomic level requires molecular dynamics simulation, which resolves the full trajectory of every atom through the cascade and recovery.

This work uses LAMMPS to perform a complete cascade simulation in bulk wurtzite GaN, including equilibration, a 5 keV Ga PKA cascade, and a dynamic recovery phase at 300 K. The results are analyzed in terms of energy evolution, atomic displacements, and residual defect counts. A detailed analysis of the thermodynamic dashboard and point defect evolution is provided to quantify the damage and self healing phases.


2. Wurtzite Gallium Nitride: Relevant Properties

2.1 Crystal Structure

The wurtzite crystal structure (space group P6_3mc) consists of alternating hexagonal planes of Ga and N atoms stacked along the c axis. Each Ga atom is tetrahedrally coordinated by four N atoms and vice versa, forming a network of strongly bonded Ga N pairs. The lattice parameters used in this simulation follow the established experimental values for GaN.

Parameter Value Source
Lattice constant a 3.189 A Experimental / DFT consensus
c/a ratio 1.626 Experimental / DFT consensus
Lattice constant c 5.185 A Derived: a x c/a
Atoms per unit cell 4 (2 Ga + 2 N) Wurtzite
Bandgap 3.4 eV (direct) Experimental
Bulk modulus ~210 GPa Experimental
Melting point ~2500 K (at pressure) Experimental

2.2 Radiation Tolerance Mechanisms

GaN displays better radiation tolerance than most common semiconductors for three interconnected reasons.

  1. High displacement threshold energies. The minimum energy needed to permanently displace a Ga atom from its lattice site is approximately 25 eV, and for an N atom it is approximately 20 eV. These values are higher than in Si (about 15 eV) and GaAs (about 10 eV for As). A higher threshold means fewer atoms are displaced per unit of radiation dose.
  2. Strong Ga N bonding. The Ga N bond has a bond dissociation energy of approximately 8.9 eV, reflecting its mixed covalent and ionic character. This strong bonding creates a large restoring force on any displaced atom, driving it back toward a lattice site and facilitating Frenkel pair recombination.
  3. Compact wurtzite geometry. The short nearest neighbor distances in wurtzite GaN mean that an interstitial atom is never far from a vacancy. Close proximity increases the probability of recombination during the thermal spike and recovery phases.

3. Simulation Methodology

3.1 Software

All simulations were performed with LAMMPS (Large scale Atomic/Massively Parallel Simulator), version 22 July 2025 Update 4. Parallelization was achieved using 4 MPI tasks for the cascade run and 4 MPI tasks for equilibration.

3.2 Simulation Cell

The simulation box was constructed by repeating the 4 atom wurtzite GaN unit cell 24 times along the a1 axis, 24 times along a2, and 16 times along a3. This gives a supercell containing 55,296 atoms in a box of approximately 113.5 x 65.6 x 82.8 angstrom (after NPT relaxation). Periodic boundary conditions were applied in all three directions, eliminating surface effects and simulating bulk behavior.

3.3 Interatomic Potential

Atomic interactions were described using a hybrid overlay of two potentials.

  • Tersoff potential (GaN.tersoff): This potential governs bonding at normal interatomic separations. It captures the equilibrium lattice constants, elastic constants, surface energies, and phonon dispersion of GaN with good accuracy. The neighbor cutoff radius is 5.1 angstrom.
  • ZBL potential (Ziegler Biersack Littmark): This potential describes the strongly repulsive nuclear nuclear interaction when two atoms approach closer than about 1.2 angstrom, which occurs during high energy collisions in the cascade. ZBL is a universal screened Coulomb potential.

3.4 Three Stage Simulation Workflow

The simulation was executed as three sequential phases.

  • Phase 1: Ballistic Cascade (0-5000 steps). A Ga atom at the center of the simulation box was selected as the PKA and given a kinetic energy of 5,000 eV. The outer shell of atoms was coupled to a Langevin thermostat at 300 K. An adaptive timestep was employed, dynamically adjusting between roughly 0.00008 ps and 0.00020 ps to safely resolve high speed collisions.
  • Phase 2: Dynamic Recovery (5000-25000 steps). After the ballistic cascade, the adaptive timestep was replaced with a fixed timestep and the entire box was coupled to an NVT thermostat at 300 K to simulate room temperature thermal recovery.
  • Phase 3: Quench and Minimization (25000-35000 steps). The structure was rapidly cooled to suppress thermal vibrations, allowing the permanent remaining defects to be counted without ambiguity.

4. Results

4.1 Point Defect Evolution (Wigner Seitz Analysis)

 

W-S ANALYSIS DEFECT COUNT
WIGNER SEITZ ANALYSIS - DEFECT COUNT PER FRAME PLOT

The evolution of point defects during the cascade and recovery phases is clearly illustrated in the provided Wigner Seitz analysis plot. The graph displays a near perfect overlap between the interstitial population (blue solid line) and the vacancy population (red dashed line). This exact overlap confirms that the primary mode of radiation damage in this simulation is purely stoichiometric Frenkel pair generation, with no anomalous loss of atoms.

Upon PKA injection, the defect count rises exponentially from 0, reaching an absolute peak of approximately 490 Frenkel pairs near frame 24. This zenith corresponds to the maximum spatial expansion of the ballistic cascade before the system has time to react. Immediately following this peak, the highly localized kinetic energy initiates a rapid dynamic recovery phase. A steep, continuous decline in the defect count is observed as the lattice naturally heals, with closely spaced Frenkel pairs spontaneously recombining. By frame 140, the system stabilizes, leaving a permanent residual population of approximately 113 defects.

4.2 Thermodynamic and Energy Dynamics (LAMMPS Dashboard)
 

The four panel LAMMPS Cascade Simulation Dashboard captures the exact physical state of the crystal across the entire 35,000 step timeline.

  • Thermal Spike and Quench (Top Left): Immediately following the 5 keV PKA strike, the localized kinetic energy creates an intense thermal spike, driving the effective system temperature to roughly 1,000 K. The Langevin thermostat boundaries successfully act as a bulk material heat sink, rapidly draining this excess energy and pulling the temperature down to the 300 K baseline just prior to Phase 2. The temperature remains completely stable during the NVT dynamic recovery window, before smoothly decaying to roughly 100 K during the final Phase 3 structural quench.
  • Energy Evolution (Top Right): The structural damage is directly mirrored in the potential energy (blue curve). It surges dramatically from its 0 K reference state (approximately −245,800 eV) up to a peak near −241,800 eV, representing the massive internal strain and bond breaking of the cascade. As the Frenkel pairs recombine and the lattice relaxes, the potential energy drops and plateaus near −243,500 eV during the room temperature NVT phase. The total energy (green curve) steadily decreases over the entire simulation as the thermostats remove the foreign 5 keV energy packet injected by the initial radiation strike.
  • Shockwave Dissipation (Bottom Left): The violent primary collision generates a massive compressive stress wave radiating outward from the impact zone. This registers as a severe pressure spike of approximately 26,500 bar near step 1000. This acoustic shockwave dissipates rapidly. By the time the system enters the Phase 2 recovery window, the pressure stabilizes at roughly 7,500 bar. This remaining elevated pressure is not a dynamic wave, but rather the static residual elastic stress field exerted by the surviving 113 permanent defects wedged into the crystal lattice.
  • Adaptive Timestep Optimization (Bottom Right): The adaptive timestep plot confirms the numerical stability of the ballistic phase. During the most chaotic, high velocity initial collisions (steps 0-600), the solver aggressively restricted the timestep down to roughly 0.00008 ps to prevent unphysical atomic overlaps. As the kinetic energy spread out and the atoms slowed down, the timestep safely scaled back up to its designated 0.00020 ps ceiling.

4.3 Quantification of Self Healing

The data extracted from the visual graphs allows for a precise quantification of the wurtzite crystal's self healing efficiency.

Recovery Metric Value
Peak Frenkel pairs (maximum damage) ~490
Surviving Frenkel pairs (permanent damage) ~113
Recombined Frenkel pairs ~377
Defect Recovery Fraction ~76.9 %
Peak Cascade Pressure ~26,500 bar
Peak Thermal Spike Temperature ~1,000 K

5. Conclusion

The simulated 5 keV PKA cascade vividly demonstrates the intrinsic radiation hardness of gallium nitride. The quantitative Wigner Seitz data reveals an exceptional instantaneous self healing capability: of the roughly 490 Frenkel pairs generated at the absolute peak of the ballistic cascade, approximately 77 percent successfully recombined and healed within mere picoseconds at room temperature, leaving only 113 permanent defects.

This profound resilience is governed by fundamental atomic physics. The strongly mixed ionic covalent nature of the Ga N bond acts as a powerful restoring spring. When combined with the highly compact nature of the wurtzite crystal lattice, displaced interstitial atoms are physically restricted from migrating far from their parent vacancies. Consequently, the brief, intense heat generated by the 1,000 K thermal spike provides the exact localized activation energy required to snap these nearby displaced atoms back into their correct crystallographic sites before the crystal cools.

These atomic level dynamics validate why GaN components consistently outperform traditional silicon or gallium arsenide semiconductors in harsh radiation environments. Because the material dynamically anneals out the vast majority of primary collision damage at standard operating temperatures, cumulative device degradation is drastically slowed. This intrinsic self healing positions gallium nitride as a premier, indispensable material for next generation spaceflight electronics, high altitude avionics, nuclear reactor diagnostics, and high energy particle accelerators where long term structural and electrical reliability is an absolute requirement.

Report Information

Explore More Projects

View All 2026 Projects