Charge-optimized many-body interaction potential for AlN revisited to explore plasma–surface interactions

Plasma–surface interactions during AlN thin film sputter deposition could be studied by means of reactive molecular dynamics (RMD) methods. This requires an interaction potential that describes all species as well as wall interactions (e.g., particle emission, damage formation) appropriately. However, previous works focused on the establishment of AlN bulk potentials. Although for the third-generation charge-optimized many-body (COMB3) potential at least a single reference surface was taken into account, surface interactions are subject to limited reliability only. The demand for a revised COMB3 AlN potential is met in two steps: First, the Ziegler–Biersack–Littmark potential is tapered and the variable charge model QTE\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{+}$$\end{document}+ is implemented to account for high-energy collisions and distant charge transport, respectively. Second, the underlying parameterization is reworked by applying a self-adaptive evolution strategy implemented in the GARFfield software. Four wurtzite, three zinc blende and three rock salt surfaces are considered. An example study on the ion bombardment induced particle emission and point defect formation reveals that the revised COMB3 AlN potential is appropriate for the accurate investigation of plasma–surface interactions by means of RMD simulations.


Introduction
The powering of micro-electro-mechanical systems and wireless sensors can be realized by energy harvesting (scavenging), overcoming the limitations of batteries (e.g., runtime, maintenance costs).Absorbed ambient vibrations can be translated to mechanical stresses of a piezoelectric material.The induced strain causes a polarization and, hence, voltage drop across the material.Hexagonal aluminium nitride (AlN) is argued to be a promising material system for such an application due its high chemical, thermal, and mechanical stability [1][2][3][4] .Its properties make it also a widely used subsystem for hard coatings and protective wear applications (e.g., transition metal aluminium nitride, transition metal aluminium oxynitride) 5 .For either application (i.e., energy harvesting, hard coatings) sputter deposition is exploited to fabricate the thin films.A plasma is ignited between two opposing electrodes by applying a voltage of up to 1 kV.Background gas atoms (e.g., Ar, N, N 2 ) are ionized and accelerated toward the (Al) target electrode.The ion bombardment spawns a collision cascade within the material which eventually leads to the emission of surface atoms.After transport through the gas phase these atoms may impinge onto the opposing substrate surface.The simultaneous flux of neutral background gas atoms and (less energetic) ion bombardment on the substrate surface together result in sputter deposited (AlN) thin film growth.The film properties (e.g., stress, defect structure) intrinsically determine the performance during application.Hence, it is fundamental to truly understand and possibly design the defect formation during deposition.
The required information can accessed by reactive molecular dynamics (RMD) simulations of AlN growth or sputtering processes.However, the fidelity of RMD studies inherently relies on the accuracy of the applied interaction potential.A variety of potentials (e.g., Vashishta, Tersoff) has been proposed for the investigation of AlN bulk systems and processes (e.g., pressure induced phase changes) 6,7 .A detailed comparison can be found elsewhere 8 .Surfaces have only been taken into account for setting up the third-generation charge-optimized many-body (COMB3) AlN potential, which is also able to describe Al and N 2 9, 10 .However, the underlying database consisits of only a single surface energy and four point defect formation energies.This shortcoming may be of little relevance for the intended application (i.e., nanostructures, heterogeneous interfaces), but the lack of relevant reference data may cause issues for growth or sputtering simulations 9 .Notably, the published COMB3 TiN, TiO 2 and Al 2 O 3 potentials may provide a basis for developing COMB3 potentials for hard coatings such as TiAlN or even

Interaction potential
In the frame of COMB3, a system that consists of two atoms of one kind, is described by short range and van der Waals interactions.The latter is often neglected due to its relatively weak contribution.When two atoms i and j differ in their elemental type, charge is transferred from on to the other.The charge of each atom consists of the particular nucleus charge Z, which is modeled as a point charge, and q − Z, which is spatially described by a 1s Slater type orbtial (STO).If required, point dipoles can be taken into account to model the system's electrostatic interaction more accurately.The charge transport takes place on the electronic time scale.As a consequence of the mass ratio, the nucleus can be thought of as immobile during the electron dynamics.In combination with COMB3, the corresponding equilibration process is typically modeled by the charge equilibration (QEq) method in conjunction with an extended Lagrangian formalism 13,16 .However, QEq fails to describe the charge transport during most relevant surface interactions (e.g., adsorption, desorption, deposition, sputtering).The application of the charge transfer equilibration (QTE + ) in favor of the QEq model is thus recommended for studies that target surface interactions 17 .The 1s STOs are required to enable a matching charge transfer per unit time.The corresponding electrostatic long range interactions are described in detail elsewhere 13,16 .
For a binary system, the initially mentioned short range interactions consist of a repulsive as well as an attractive pair function V A and V R , respectively 13,16 : r i j is the distance between atom i and atom j.A detailed description of the per element parameter λ * i , D i (q i ), D j (q j ), α * i and B i j * (q i , q j ) can be found in 13,16 .Only the pair parameters A i j , λ i j , B 0 i j and α 0 i j are modified for Al-N interactions in this work (B 1 i j = B2 i j = α 1 i j = α 2 i j = 0).The attractive term V A is scaled with the average bond order of both atoms b i j+b ji 2 which allows for a more flexible and environment dependent potential 13,16 .The bond order b i j in a carbon exclusive system (i.e., AlN) is defined as follows f c is the Tersoff cutoff function.ζ i j = exp(ω i j (r i j − r ik )) penalizes bond asymmetries.The parameter ω i j is used to scale its impact 13,16 .It is adjusted in this work for N-N, Al-N and N-Al pairs.The coefficients b n,i jk of a sixth order polynomial are used to introduce a dependency on the bond angle θ i jk .The i-th atom is the centering atom.These parameters are three dimensional, whereas in the original COMB3 formalism they are intended to be two dimensional (i.e., b n,i j ) 13,16 .They are fitted in this work for Al-Al-N, Al-N-Al, Al-N-N, N-Al-Al, N-Al-N and N-N-Al.The effect of the coordination number Ω i is added by P i j (Ω i ) = c 0,i j Ω i + c 1,i j exp(c 2,i j Ω i ) + c 3,i j for Al-N and N-Al.
Legendre polynomials up to the sixth order are added as a function of cos(θ i jk ) to the potential energy to ease the differentiation between face-centred cubic (fcc) and hexagonal close-packed (hcp) phases 13,16 .The corresponding coefficients i jk are adjusted for Al-Al-N (Al-N-Al), Al-N-N, N-Al-Al and N-Al-N (N-N-Al).A symmetry is imposed that demands K LP1 i jk = K LP1 ik j .High-energy collisions during the ion bombardment are accounted for by implementing a combined (tapered) version of the revised COMB3 2022 QTE + with the Ziegler-Biersack-Littmark (ZBL) potential in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [18][19][20][21] .For r th,2 < r i j , the interaction is purely described by the COMB3 2022 QTE + potential.In between r th,1 < r i j ≤ r th,2 , the interactions are determined by the tapered interaction potential.The contribution of the ZBL potential is gradually increased with decreasing distance by scaling with the Tersoff cutoff function.For r i j ≤ r th,1 , the interaction is purely described by the ZBL potential.The values for r th,1 and r th,2 for the considered element pairs are listed in Table 1.Bond-order definitions and calculations are not changed by the implemented tapering.Hence, two atoms with r th,1 < r i j ≤ r th,2 contribute to each other's bond-order.A similar tapering of the Lennard-Jones (LJ) potential (with an outer cutoff radius of 8.5 Å) with the ZBL potential is used for Ar-Ar interactions 18 .The values for r th,1 and r th,2 are also listed in Table 1.The depths of the LJ potential well and zero-crossing distance are 0.0104 and 3.4 Å, respectively.The Al-Ar and N-Ar interactions are described purely by the ZBL potential (with an outer cutoff radius of 8.0 Å).

Al-Al
The implementations are verified by comparing the integrated force with the corresponding potential energy and the tapered with the non-tapered potentials for a relevant set of cases (see Supplementary Methods).

Reference data
The assembled reference data set used to revise the COMB3 2016 * AlN potential is listed in the following: • 28 Al x N y clusters with up to eight atoms (i.e., x ≤ 6, y ≤ 4) (resembling structures on top a surface or during ion bombardment induced surface reconstructions).They are described by the atoms' geometry (i.e., bond lengths, bond angles), and (relative) binding energies.
• Five (seven) point defect structures for the wurtzite (zinc blende) phase.They are described by defect formation energies (N rich conditions are assumed).
• Four wurtzite, three zinc blende, and three rock salt surfaces are taken into account.They are described by their surface energies.
• 150 bond lengths and angles.They describe the atom geometries introduced above.

GARFfield extensions
The software GARFfield developed by Jaramillo Botero et al. is a reactive force field optimizer that is described in detail in 14 .
In the following a brief summary focusing on implemented extensions, minor revisions, and the herein applied features is presented.
Multiple interaction potentials (e.g., ReaxFF, eFF, COMB) are included in GARFfield.The COMB potential and the QEq method are replaced in this work with the COMB3 potential and the QTE + method, respectively.The reading (writing) of the interaction potential files are revised accordingly.The COMB3 AlN parameters (genes) to be fitted are listed in the preceding section.This list is referred to as a chromosome with length n in the frame of a genetic algorithm (GA) as it is applied within GARFfield 14 .An initial population with diverse chromosomes and genes is established by sampling randomly in a given interval around the COMB3 2016 * QEq AlN parameters.The lower and upper bound of the i-th gene are x lo,i and x hi,i , respectively.The fitness of each chromosome is evaluated by performing molecular statics calculations of the atom configurations mentioned in the preceding section, with the LAMMPS library (including COMB3 QTE + ) 19,20 .A conjugate gradient descent algorithm is used to relax the individual structures.The time step for the QTE + fictitious charge motions is chosen as 10 −2 fs.The tolerance for the maximum charge force and force on any atom is 10 −2 V and 1 eV/Å, respectively.These evaluations are handled in a parallel manner by utilizing a modified version of PGAPack in GARFfield 14,22 .A master-slave architecture is established, in which each processor performs the molecular statics calculations in serial for each individual out of the population.
The list of available reference properties (e.g., cell parameters, formation energies) is enlarged by the implementation of elastic constant and surface energy calculations.Moreover, the order of the program execution is altered to enable more consistent evaluations.First, the relaxed cell parameters are determined.Second, these cell parameters are then used to scale other atom configurations (e.g., point defect structures).
Deviation between the fitted and the reference properties are quantified by the weighted mean percentage error (WMPE).The sectional weights for the cell parameters, atom geometries, elastic constants, and energies are 100, 3, 1, 12, respectively.The weights for each entry in the cell parameters, atom geometries, and elastic constants are chosen to even out the impact of each atom structure.The weights for the energies are chosen to even out the impact of i) wurzite AlN heat of formation; ii) relative zinc blende, rock salt and CsCl heats of formation (heats of reaction); iii) binding energies of AlN clusters; iv) wurtzite point defect formation energies; v) zinc blende point defect formation energies; vi) surface energies.However, the WMPE is still ill defined for a small number of reference samples.This issue is resolved by defining a limit for the denominator value of 1 eV, 0.1 GPa and 0.1 eV/Å 2 for binding (formation) energies, elastic constants, and surface energies, respectively.
As a metric for the fitness of the evolution, the WMPE is used to rank and select parents for reproduction and, hence, creation of the next population.The new chromosomes are obtained by applying two-point-crossover recombination to the parents' chromosomes.Genes are then randomly selected per chromosome for mutation.The probability is commonly chosen as the inverse of the chromosome length.Thus, mutating only one gene is the most probable event 23 .The value x i of the i-th gene is updated by x i → σ N i (0, 1)x i , whereas σ and N i (0, 1) are the isotropic mutation strength and the standard normal distribution, respectively.
On the one hand, two-point-crossover recombination is found to hinder substantial evolutionary progress.This is reasoned by the specific linkage and setup of the COMB3 potential parameters.On the other hand, a GA without crossover recombination is considered to be just an inefficient evolution strategy (ES).Hence, we implemented an anisotropic self-adaptive evolution strategy with intermediate recombination (µ/µ I /λ )-σ SA-ES in GARFfield 14 .µ and λ are the number of parents and population size, respectively.Anisotropic self-adaptation means that the formerly constant mutation strength (step size) σ is replaced with an evolving mutation strength (step size) σ i for each gene i.The list of these individual mutation strengths (step sizes) is appended to the chromosome.Hence, its length n is doubled.Intermediate recombination means that the offspring population's chromosomes are at first set equal to the centroid of the parents' chromosomes.Intermediate recombination is implemented in the modified PGAPack software 22 .All genes are then mutated.First, each mutation strength σ i is updated following σ i → σ i exp(τ N (0, 1) + τN i (0, 1)).Random realizations are sampled once from τ N (0, 1) and n times from τN i (0, 1) per generation to provide the flexibility for changes on the global and individual scale, respectively.The factors τ and τ are defined by τ = 2 √ n −1 and τ = √ 2n −1 , respectively [23][24][25] .Second, each gene x i is updated by x i → σ i N i (0, 1)(x hi,i − x lo,i )/2 + x i .Boundaries for each gene demanded by the COMB3 potential formalism are enforced (e.g., no square root of negative numbers).The lower (upper) bound of the initialization interval for the i-th gene [x lo,i ,x hi,i ] is chosen to be its initial value minus (plus) 20 % when the gene directly effects either the bond-order or Legendre polynomial computations.When the gene effects only the attractive or repulsive pair potential, this percentage value is lowered by a factor of 10 (i.e., ± 2 %).
Overall the modified GARFfield implementation allows for a consistent and continuous transition from a GA to an ES, making it a more diverse evolutionary algorithm (EA).
The (1,79)-σ SA-ES implemented in GARFfield is applied to revise the COMB3 2016 * AlN potential.Ten iterations each consisting of up to 200 generations have been performed.For the final run, the range of the initialization interval is reduced by a factor of four (i.e., ± 5 %, ± 0.5 %).

Exemplary study setup: Ion impingement on AlN(0001)
As an exemplary study, the isolated impingement of ions -in the present case Al + , N + , Ar + , and N + 2 -onto the ideal wurtzite Al-AlN(0001) surface is investigated.The bombardment-induced mean damage formation and the flux emission are considered.Ion energies E ion range from 10 eV to 100 eV in 10 eV steps.This range is chosen on the basis of experimentally measured mean ion energies impinging onto the substrate during the reactive radio frequency sputtering of AlN 26 .For each ion species and energy, the simulation is repeated 50 times, whereas the ions' lateral position is sampled from a uniform distribution and scaled with the surface dimensions.The damage formation is quantified by means of the induced point defect populations (i.e., LAMMPS is used for all simulations outlined in this subsection 19,20 .First, a bulk 7 × 4 × 7 wurtzite AlN system with lattice constants a = 3.154 and c = 1.6a is constructed.It is shifted and cleaved in [001] direction to obtain the Al terminated 4/14 QEq are compared with DFT and experimental findings [28][29][30][31][32][33][34][35] .The references are broken down in the Supplementary Table S1.Labels in brackets allow for a structure identification in the provided reference when necessary.The subtraction of two labels indicates a relative binding energy. Al-AlN(0001) surface.Bottom layers of the surface slab are excluded from time integration up to a height z th to anchor it in the simulation domain, reduce the effect of the N-AlN(0001) surface on the damage formation and reduce the computational cost.
The variable z th is defined through E ion as follows z th = 4c − 3c( √ E ion − √ 10 eV)/( √ 100 eV − √ 10 eV).Four (E ion =10 eV) to one (E ion =100 eV) AlN layers are kept immobile.The atom charge distribution is calculated on each time step by applying the QTE + method again with a time step of 10 −2 fs and a tolerance for the maximum charge force of 10 −2 V.The fitted STOs ζ ov for Al and N are 0.668 Å −1 and 1.239 Å −1 , respectively 17 .
This exemplary case study is performed twice to assess the necessity for equilibrating the surface prior to a particle impact: a) The mobile surface atoms are equilibrated at 300 K by applying a Berendensen thermostat to two AlN layers with z th < z < z th + 2c 27 .A damping constant and time step of 100 fs and 0.25 fs is used, respectively.The simulations are run for 100 ps.b) To obtain a surface relaxed initially at 0 K, atom positions are adjusted by applying a conjugate gradient minimization until the force on any atom is below 1 eV/Å.For the latter case, a lower bound for the final temperature can be predicted N and k B are the number of active (mobile) atoms and the Boltzmann constant, respectively.In either case the equilibrated ideal Al-AlN(0001) surfaces are then reloaded for each ion bombardment.A consecutive bombardment inherently leads to a drifting surface state.A quantification and analysis of the overlaying effects of different ion energies, ion species, and surface states is out of scope for this work, however.
To simulate the ion bombardment, the time step is eventually lowered starting at 0.25 fs to secure a maximum change of the atoms' kinetic energy and displacement below 0.1 eV and 0.01 Å, respectively.The simulation time t per impingement is defined by t = 2 ps + 3( √ E ion − √ 10 eV)/( √ 100 eV − √ 10 eV) ps, to reflect the longer relaxation time for higher ion energies.2 ps (E ion =10 eV) to 5 ps (E ion =100 eV) long RMD simulations are performed for each energy.

Results
First, the results of the COMB3 2016 * QEq and the revised COMB3 2022 QTE + AlN potential on the reference set outlined previously are compared to density functional theory (DFT) and experimental (exp.)findings.Second, an examplary study of the ion bombardment induced flux emission and the damage formation is presented.

Performance on reference set
The fidelity of modeling Al x N y cluster is assessed at first.The (relative) binding energies for clusters with up to eight atoms (x ≤ 6, y ≤ 4) are shown in Fig. 1.Dissociation energies are converted to binding energies by imposing a binding energy of -5.06 eV/atom and -0.77 eV/atom for N 2 and Al 2 , respectively 29 .COMB3 2016 * QEq describes the differentiation between exothermic and endothermic binding energies correctly.However, the absolute values are found to deviate significantly from the reference values.In particular the binding energy of Al 1 N 1 are approximately three times too high.This could introduce exaggerated sputtering of such monomers due to ion bombardment.An approximately threefold overestimation is also found for Al 1 N 2 (N-Al-N,D ∞h 31 ).In general it is found that COMB3  Furthermore, it struggles with the v Al and the (N-N) N point defect structures presented in Fig. 2. The negative values are likely to cause an exaggerated number of vacant Al sites and incorporated N split interstitials during an ongoing ion bombardment.In contrast, COMB3 2022 QTE + is found to describe this and the other point defects reasonably well.Both potentials describe the order of surface energies shown in Fig. 3 S2 and Supplementary Table S3.constant C 44 accurately.Though both models also under-and overestimate C 12 and C 11 , respectively.The reference heat of formation ∆H f and heat of reaction ∆H rxn match reasonably well.
Concerning the point defect structure shown in Fig. 2, COMB3 2016 * QEq describes only Al i and N Al accurately.For all other point defects, deviations of more than 2 eV are observed.The binding energy of vacant Al sites v Al and N split interstitials (N-N) N even deviate by more than 178 % and 177 %, respectively.This issue with these two point defect structures is also observed for the wurtzite structure, discussed in the preceding paragraph (see Table 2).This repetition signifies a systematic problem of the COMB3 2016 * QEq potential with v Al and (N-N) N .The last is assumed to be caused by its nitrogen bond configuration resembling molecular nitrogen.In contrast, COMB3 2022 QTE + describes all point defect structures with high fidelity.Both models describe the order of surface energies shown in Fig. 3    The properties considered for the rock salt structure are listed in Table 4. COMB3 2016 * QEq overestimates the lattice constants and fails to resolve the heat of formation and the heat of reaction.In contrast, COMB3 2022 QTE + describes the lattice constant, the heat of formation, and the heat of reaction accurately.
However, both models are not capable of describing the correct order of surface energies depicted in Fig. 3, whereas COMB3 2022 QTE + provides a smaller mean absolute error (i.e., 0.051 eV/Å 2 ) than COMB3 2016 * QEq (i.e., 0.096 eV/Å 2 ).The heat of reaction of CsCl relative to the wurtzite structure ∆H CsCl−WZ rxn for COMB3 2022 QTE + , COMB3 2016 * QEq , COMB3 2016 QEq , and DFT reference are 4.02 eV, 1.37 eV, 1.37 eV, and 4.03 eV, respectively 9 .Hence, COMB3 2022 QTE + provides a more accurate description of the phase stability.Relaxed atom structures resemble the referenced configurations for all herein considered cases by means of bond lengths and bond angles (not shown).QEq AlN publication, and DFT as well as experimental findings 9,39,47 .The references are broken down in the Supplementary Table S2, Supplementary Table S3, and Supplementary Table S4.

Exemplary study: Ion impingement on AlN(0001)
In the following, the irradiation of an ideal wurtzite AlN(0001) surface by Al + , N + , Ar + , and N + 2 ions is investigated.The emitted nitrogen flux Γ out N per incident ion flux Γ in ion is displayed in Fig. 4 (a).As evident, a reflection probability for N + of 20.0 % for E ion = 10 eV decreases down to 12.0 % for E ion ≥ 60 eV (root-mean-square deviation: 2.2 %).Above 20 eV a conversion of incident nitrogen ions and surface atoms to molecular nitrogen is observed.This process peaks at 40 eV as depicted in Fig. 4  (b).An incident nitrogen ion eventually hits a surface nitrogen atom and forms a nitrogen split interstitial (N-N) N at the surface for a brief moment until the inverted momentum (now in surface normal direction) leads to the emission of the just formed (N-N) N as a N 2 molecule.For even higher ion energies (i.e., E ion > 40 eV) this process becomes less relevant due the more likely and deeper ion implantation as depicted in Fig. 5.
The same process also leads to an increase in the emission of molecular nitrogen for incident N + 2 that is split up into two nitrogen atoms at the surface for ion energies of 40-80 eV.This process is peaked at 60-70 eV.The broadening and shift to higher ion energies is reasoned by a decreased implantation depth (see Fig. 5).This is due the initially shared and eventually separated momentum.The implantation depths of N + and N + 2 for approximately 40 eV and 60-70 eV, respectively, are indeed comparable.Hence, the impingement of N + 2 intrinsically affects a region closer to the surface which facilitates the formation and eventual emission of split interstitials as molecular nitrogen.
The incoming flux of N + 2 contributes most significantly to the outgoing nitrogen atom flux for E ion > 20 eV due to a splitting and subsequent reflection of one N atom at the surface.This process increases only slightly for energies E ion > 40 eV.The remaining second N atom may be implanted, stick, reflected, or leave the surface with another N surface atom as an N 2 molecule.The last process contributes most substantially to N 2 emission for energies of 30-80 eV (sticking is dominant at lower energies and implantation is dominant at higher energies).
Only a negligible number of Al ions and Al surface atoms are reflected and sputtered, respectively (not shown).This is due to high sticking probabilities and small sputtering yields in the considered ion energy range (i.e., 10-100 eV).The incorporated Al interstitials Al i are quantified by means of their defect population as shown in Fig. 6 (b).Up to approximately 30 eV, solely the implantation of impinging Al + ions contributes to the Al interstitial population.For higher ion energies, forward sputtering (peening) of Al surface atoms leads to a steadily increasing number of Al interstitials for any in this work considered ion species.
All Ar ions are reflected for ion energies equal to 10-20 eV.At 30 eV 76 % of the incident Ar ions are still reflected.The reflection probability decreases linearly with increasing ion energies down to 29 % at E ion = 100 eV.The N and Al vacancy populations are increased for higher ion energies due to the forward sputtering (peening) of N and Al atoms, respectively (see  2 ).Solid and dashed lines indicate surfaces initially equilibrated at 300 K and relaxed at 0 K, respectively.Supplementary Fig. S6).Only a negligible amount of other point defect populations (i.e., ρ Al N , ρ N Al , ρ Ar Al , ρ Ar N ) are observed (not shown).
Starting off with a relaxed surface at 0 K (see equation ( 3)) provides comparable results for nearly all properties (e.g., Ar reflection probability, implantation depth).However, a noticeable deviation is observed for the emitted molecular nitrogen flux Γ out N 2 as a function of the incident N + 2 ion flux Γ in ion for 40 eV< E ion < 80 eV, as displayed in Fig. 4 (b).The described process of splitting, forming split interstitials, and leaving the system as N 2 molecules is weakened when starting with a relaxed surface at 0 K.A too large proportion of the ions' kinetic energy is transferred to the surface atoms which eventually hinders this process for this ion energy range.For higher ion energies (i.e., 90-100 eV) this issue is resolved due to the deeper implantation depth (see Fig. 4 (b) and Fig. 5).Moreover, it is found that starting with a relaxed surface at 0 K promotes the ion implantation and generation of corresponding interstitials (see Al + in Fig. 6 (b), N + in Fig. 6 (a), and N + 2 in Fig. 6 (a)).It has a minor effect on the ion bombardment induced Frenkel pair generation.
The corresponding surface slab temperature T is depicted in Fig. 7. T sim and T pred are the simulated mean temperature and its lower bound as predicted by theory (see equation ( 3)).A good agreement between theory and simulation is observed.The first may be used to design studies starting with a relaxed surface at 0 K and, hence, bypassing the burden of equilibrating a surface at a specific temperature.However, the provided theoretical estimate should be used with care when the initial atom configuration is high in potential energy (relative to its ground state).

Conclusions
GARFfield with its new extensions allowed for the revision of the COMB3 2016 QEq AlN potential.COMB3 2022 QTE + is proposed to enable RMD studies of the wurtzite and zinc blende phase AlN plasma-surface interactions.The potential overcomes manifold shortcomings for the intended application (e.g., negative defect formation energies for v Al and (N-N) N ) of the COMB3 2016 QEq potential.The excellent agreement of the results obtained using the COMB3 2022 QTE + potential with experimental references and DFT based computations indicates that the damage generation (e.g., Frenkel pair) during an ongoing ion bombardment is appropriately described.High-energy collisions are taken into account by merging (tapering) the COMB3 2022 QTE + with the ZBL potential.The decent descriptions of most Al x N y clusters with up to eight atoms (i.e., x ≤ 6, y ≤ 4) demonstrate that the COMB3 2022 QTE + potential may be considered for systems involving under-coordinated atoms (e.g., atop surfaces).Moreover, it is expected to describe more reliably the relative amounts of different sputtered species (e.g., N, N 2 , Al, Al 1 N 1 , Al 2 ).
The COMB3 2022 QTE + may, however, be used with care when studying AlN in the rock salt phase.The potential is found to describe the lattice parameter and heat of formation (reaction) with high accuracy.But the correct order of surface energies could not be reproduced and, most importantly, only a small set reference properties is considered for this phase.
An example study on the effect of the isolated ion bombardment (i.e., Al + , N + , Ar + , N + 2 ) onto the ideal wurtzite Al-AlN(0001) surface was conducted.An analysis of the individual trajectories, particle emissions, and point defect formations for 10/14 singular impingement events revealed that N + and N + 2 may form N split interstitials with surface N atoms.These eventually leave as N 2 for ion energies in the range of 30-50 eV and 40-80 eV, respectively.The broadening and shift to higher energies for N + 2 is explained by the initially shared and finally separated momentum.This process is weakened when the surfaces are initially not equilibrated at the 300 K, but relaxed at 0 K.In addition, in this case the incorporation of the incident ions is facilitated.Notably, most other properties (e.g., Frenkel pair populations, implantation depths) are found to be comparable.It is argued that such a procedure may be of interest for applications where its shortcomings are less critical, but the reduced computational costs are of importance (e.g., for parameters screenings).
Self-adaptive evolution strategies are implemented into the GARFfield framework and allow for a more versatile interaction potential optimization.A continuous transition or alternation between genetic algorithm and evolution strategy can be readily applied, providing a more versatile evolutionary algorithm in GARFfield.Further extensions such as surface energy calculations or replacing the COMB with the COMB3 potential enhance the application possibilities.

Figure 4 .
Figure 4. (a) Nitrogen flux Γ out N and (b) molecular nitrogen flux Γ out N 2 per incident ion flux Γ in ion are presented as a function of the ion energy E and ion species (i.e., Al + , N + , Ar + , N + 2 ).Solid and dashed lines indicate surface initially equilibrated at 300 K and relaxed at 0 K, respectively.

Figure 5 .
Figure 5. Implantation depth is presented as a function of the ion energy E and ion species (i.e., Al + , N + , Ar + , N + 2 ).Solid and dashed lines indicate surfaces initially equilibrated at 300 K and relaxed at 0 K, respectively.

Figure 6 .
Figure 6.(a) Nitrogen split interstitial population ρ N sp i and (b) Aluminum interstitial population ρ Al i are presented as a function of the ion energy E and ion species (i.e., Al + , N + , Ar + , N + 2 ).Solid and dashed lines indicate surfaces initially equilibrated at 300 K and relaxed at 0 K, respectively.

Figure 7 .
Figure 7. Theoretically predicted (see equation (3)) and by simulation obtained surface slab temperatures T pred and T sim are presented as function of the ion energy E.
3216 * QEq describes most clusters with up to three Al and N atoms to be too energetically favorable.In constrast, COMB32022QTE + models almost all clusters appropriately.The systematic overestimation is overcome.A larger deviation of ≥ 1 eV is only found for Al 4 N 4 (Xs32).

Table 2 .
Wurtzite AlN properties obtained with COMB3 2022 QTE + and COMB3 2016 * QEq are compared to the COMB3 2016 QEq AlN publication, and DFT as well as experimental findings.The considered wutzite AlN properties are listed in Table 2.The results of COMB3 2022 QTE + and COMB3 2016 * QEq are compared with DFT and experimental findings.Both potentials describe the lattice constants appropriately, but tend to overestimate the elastic constants C 11 , C 33 and C 44 .This bias is more pronounced for the AlN potential revised in this work.The heat of formation is modeled with high accuracy with COMB3 2022 QTE + , whereas the COMB3 2016 *QEq potential underestimates this parameter by approximately 0.2 eV.

Table 3 .
Zinc blende AlN properties obtained with COMB3 2022 QTE + and COMB3 2016 * QEq are compared to the COMB3 2016 QEq AlN publication, and DFT as well as experimental findings.The considered properties for the zinc blende AlN structure are listed in Table 3.The results of COMB3 2022 QTE + and COMB3 2016 *QEq are again compared to DFT and experimental findings.Both models describe the lattice constant and elastic

Table S1 .
Al x N y cluster structure COMB32022Binding energies E bind (eV/atom) of Al x N y cluster (x ≤ 6, y ≤ 4) computed with COMB3 2022 QTE + and COMB3 2016 * QEq are compared with DFT and experimental findings.Labels in brackets allow for a structure identification in the provided reference when necessary.The subtraction of two labels indicates a relative binding energy.

Table S2 .
Wurtzite AlN defect formation and surface energies obtained with COMB3 2022 QTE + and COMB3 2016 * QEq are compared to the COMB3 2016 QEq AlN publication, and DFT as well as experimental findings.