Etch Plasma–Surface Molecular Dynamics (MD) Modeling resolves the femtosecond-to-nanosecond atomic trajectories that turn an incident ion, radical, neutral, or reaction product into reflection, adsorption, bond rearrangement, implantation, damage, sputtering, etch products, and lattice heat. Its most useful reactor-scale output is not an attractive trajectory movie. It is a validated, state-conditioned statistical closure—yields, probabilities, product multiplicities, energy-angle kernels, retained composition, and uncertainty—that feature Monte Carlo, surface kMC, and continuum profile models can consume without violating atoms, energy, charge, or probability.

This upgraded page owns atomistic impact dynamics at a material surface. Electronic-structure calculations own reference energies, forces, charge states, barriers, and reaction paths; MD owns the short-time many-atom response; surface kMC owns rare thermal events and long waiting times; feature Monte Carlo owns transport between surface encounters; level set or cellular models own micron-scale geometry evolution. A trustworthy workflow makes each handoff explicit and never asks an MD cell at an artificial flux to impersonate an entire etch chamber.

Begin with the equation and its approximation. Classical MD integrates

where the potential energy $U$ is the central model. Velocity Verlet, for example, advances

and completes the velocity after evaluating new forces. Newton's equation is not the major uncertainty; the representation of $U$, electronic effects, boundary conditions, and sampled surface ensemble usually are.

Born–Oppenheimer or ab initio MD evaluates forces with electronic structure at each step. It offers a chemically grounded reference but is limited to relatively small cells and short trajectories. Classical reactive MD uses a fitted analytical force field and reaches larger cells/ensembles, but inherits the form and training limitations of that force field. Machine-learned interatomic potentials can approach reference-force accuracy at much lower cost, yet require careful coverage, uncertainty detection, and stable behavior outside familiar configurations. Nonreactive pair or bond-order potentials may be excellent for collision cascades in a fixed chemistry and unusable for fluorination, oxidation, polymer formation, or volatile-product creation.

No label—ReaxFF, COMB, Tersoff, neural equivariant potential, Gaussian approximation potential—guarantees validity. State the exact parameterization, code/version, long-range treatment, charge equilibration settings, taper/cutoff, element combinations, and known exclusions. Two parameter sets with the same family name can predict different surface products and thresholds.

Build the potential around the decisions the model must support. A plasma-etch training/qualification set should cover more than equilibrium bulk structure. Include representative bulk and surface polymorphs; clean and terminated facets; amorphous environments; under- and over-coordinated atoms; adsorbates; radicals; molecules and ions where the model claims them; bond stretching/breaking; reaction paths; collision-compressed configurations; defects; interfaces; and volatile products.

Reference targets may include energies, forces, stresses, charges, surface/adsorption energies, reaction barriers, dissociation curves, equations of state, and product geometries. Weighting matters: an excellent bulk fit can still miss low-coordination etch chemistry. Split train, validation, and challenge sets by configuration family, not near-duplicate frames.

For configuration $c$, a generic fitting objective is

Report errors by chemical/process regime, not only one aggregate RMSE. A small average force error dominated by equilibrium atoms says little about a rare impact configuration that decides whether SiF$_4$ escapes. Test qualitative chemistry: correct stable products, bond order, barrier ordering, repulsive wall, dissociation, coordination, and heat of formation.

Use committee disagreement, descriptor distance, extrapolation grade, or a calibrated error model to flag frames for new reference calculations. Test the trigger against held-out errors. Stop or quarantine extrapolative trajectories; accepting them silently poisons the kernel.

Construct an ensemble of physical surfaces, not one ideal snapshot. Define crystal orientation and reconstruction, or generate independently quenched amorphous cells with converged density, composition, coordination and residual stress. A real plasma-facing surface may contain native oxide, halogen coverage, hydrogen, carbonaceous/polymer film, vacancies, implanted atoms, dangling bonds, mixed layers, roughness, and subsurface damage. These variables can dominate impact outcome.

Equilibrate the specimen at the intended surface temperature using a documented ensemble. Confirm energy, temperature, stress, density, radial/angle distributions, coordination and surface structure. For amorphous materials use multiple independent melt-quench or experimentally constrained replicas; frames taken a few timesteps apart from one trajectory are not independent material samples.

Use periodic boundaries laterally only when the feature represented is compatible. The lateral box must exceed the cascade/product correlation length and the potential cutoff with margin; an ejecta or pressure wave must not interact with its periodic image. The slab must be thick enough that collision/damage and heat propagation do not reach the constrained bottom before the outcome window closes.

A common slab partition has a bottom fixed/support layer, a heat-removal layer, and an upper Newtonian impact region. Constrain only enough to prevent center-of-mass drift and represent bulk support. Thermostatting atoms in the collision zone removes the very energy and reaction dynamics being measured. Converge all layer thicknesses and show the cascade is separated from the boundary.

Surface preparation is part of the model. Specify adsorption sites/correlations, polymer composition/density/cross-linking, and oxide stoichiometry/interface. Relax without erasing intended nonequilibrium state, then archive coordinates, velocities, charges, history, and seed.

Sample incident particles from the upstream distribution. An impact ensemble is conditioned on incident species $s$, kinetic energy $E$, direction $\Omega$, surface state $\chi$, material $m$, and temperature $T_s$. If the downstream closure will be used over a range, sample a designed grid or weighted distribution across that range.

For an incident particle of mass $M$ and kinetic energy $E$,

Transform direction from the feature or wafer frame into the local surface-normal frame. Declare whether polar angle is measured from the normal and whether angular density is per degree, per polar angle, or per solid angle. Uniform polar angle is not isotropic solid angle. If a sheath model supplies a correlated IEDF–IADF $p(E,\Omega\mid s,t)$, preserve that correlation when constructing impacts.

Randomize lateral impact position over the irreducible surface cell, accounting for symmetry only when the surface state shares it. Randomize molecular orientation and, where relevant, rotational/vibrational state. Radical and ion electronic states are difficult in ground-state classical MD; state the approximation rather than encoding an unsupported label.

Launch above the interaction range with no initial overlap. Verify the initial potential energy and force are negligible or physically intended. A launch that is too high wastes time; too low injects arbitrary potential energy. For long-range Coulomb models, "outside the cutoff" may not exist, so define reference energy and electrostatic boundary consistently.

Charged projectiles expose a fundamental boundary of ordinary MD. Fixed-topology classical MD generally does not represent neutralization, electron emission, electronic excitation, dynamic charge transfer, image charge, or plasma sheath fields faithfully. Charge-equilibration models approximate some redistribution but are not a universal electron dynamics solution. Ab initio nonadiabatic/electronic approaches may be needed for a narrowly scoped question.

If the ion is initialized as a charged atom and the potential lacks a validated charge mechanism, do not claim charge-state-resolved chemistry. It may be more defensible to treat the atom as neutral after near-surface neutralization while taking incident kinetic distributions from the sheath, then separately account for charge/current. Document the dividing surface and compare alternatives.

At higher energy, electronic stopping can remove energy not transferred to nuclei. A qualified friction model may use

with removed work tallied electronically. Do not double count it or apply it indiscriminately to thermal atoms. Nuclear stopping also requires a validated short-range repulsive wall and smooth joining.

Choose the timestep for the fastest collision, not the equilibrium vibration. A 1 fs step that conserves energy in a room-temperature crystal may fail during a 500 eV close encounter. Use a smaller fixed step or a carefully verified adaptive scheme based on maximum displacement, acceleration, or energy drift. Converge yield/product classifications with timestep, not merely average temperature.

Run a pre-impact NVE interval and monitor total energy/momentum. During impact, maintain a ledger of incident kinetic/internal energy, potential change, outgoing product energy, lattice kinetic/strain energy, thermostat work, electronic stopping, and numerical residual. The resolved ledger should satisfy

under the model's energy convention. Reaction energy is already part of a consistent potential surface; do not add it twice.

Thermostats represent heat flow to material outside the finite cell. Apply them sufficiently far from the active surface after—or weakly during—impact only when convergence supports it. Compare NVE impact windows with Langevin, Nosé–Hoover, and/or boundary damping choices. Report thermostat work. Strong global velocity rescaling can quench products, suppress sputter, and fake temperature control.

Fixed-layer momentum reflection and returning phonons alter late outcomes. Converge thickness/time, use qualified absorbing layers, or finish analysis before echoes return.

Define events from atom identities and persistence. A trajectory must end long enough after impact to distinguish transient displacement from product escape, while avoiding boundary echoes and artificial recombination. Event rules need spatial surfaces, bonding/cluster criteria, direction, energy, and persistence time.

Adsorption means the incident species remains bound in/at the surface under a defined residence observation—not simply that it has not escaped at 2 ps. Reflection means it exits the analysis plane without the intended reaction, with species/cluster, energy and direction recorded. Implantation means retained depth beyond a declared surface reference. Damage includes stable vacancies, interstitials, coordination changes, mixing or amorphization after thermal relaxation.

Physical sputter yield counts substrate atoms or formula units ejected primarily by momentum transfer. Chemical etch yield counts material removed in reaction products and must define the material unit. Total removal may combine them, but retain causal channels. A product crossing the escape plane and later returning under periodic/electrostatic conditions should not be double counted.

Use chemistry-compatible cluster analysis. Because distance bonds flicker and bond-order thresholds are model-specific, apply persistence and test against inspected trajectories. Define mixed and nonvolatile fragments.

For incident history $p$, keep a per-event atom ledger

element by element. Periodic atoms, fixed reservoirs, deposited atoms and deleted escaped particles remain traceable by immutable IDs. Report ledger failures rather than repairing them in post-processing.

Estimate yields as random variables. If $n_p^{(j)}$ is the multiplicity of outcome/product $j$ in history $p$ with statistical weight $w_p$,

Reflection or adsorption probability is the same estimator with a binary indicator. Yield can exceed one; it is not clipped into $[0,1]$. Track sample variance, covariance among products, effective sample size, and confidence intervals appropriate to binary, count, or heavy-tailed outcomes.

Independent replicates should span impact location, thermal velocities, amorphous/surface configurations, molecular orientation, and stochastic potential settings where applicable. Multiple impacts on one evolving slab are physically dose-correlated and cannot be treated as independent Bernoulli shots. Use hierarchical/cluster bootstrap or independent bombardment replicas.

Near thresholds, zero observed events does not prove zero probability. For $N$ independent unweighted trials with zero events, the approximate 95% upper bound is about $3/N$. Adaptive sampling can focus near threshold or rare tails, but likelihood weights and stopping rules must be recorded.

Resolve convergence along separate axes: lateral/vertical cell size, initial surface ensemble, timestep, trajectory duration, thermostat/boundary, impact points, energy/angle bins, number of replicas, classifier threshold, and dose. Thousands of trajectories on one undersized slab reduce sampling noise while preserving model bias.

Export full conditional kernels when transport needs correlations. A feature model may need the outgoing product/reflection measure

where $\mu$ is multiplicity or associated retained-state increment. Its integral is an outcome probability or expected multiplicity, depending on definition. State the measure $dE'\,d\Omega'$ and normalize every conditional table.

Energy and angle are often correlated: a grazing projectile may retain energy and reflect forward, while a penetrating event deposits energy and produces broad low-energy fragments. Exporting separate energy and angle marginals destroys this coupling. Species/product multiplicity may also correlate with retained damage and surface composition; preserve joint samples or a generative/sparse representation where downstream sensitivity requires it.

For outcome probability kernels,

For multiplicity kernels the integral gives $Y_j$, so the sum may exceed one. Mark which convention applies. Provide bin edges, coordinate frame, local-normal definition, interpolation, extrapolation rule, support, uncertainty, and provenance hash.

Interpolate probabilities in a positivity- and normalization-preserving parameterization such as logits/simplex methods. Interpolate yields nonnegatively and avoid smoothing away threshold behavior without evidence. Never extrapolate reactive potentials or kernels far beyond trained energy, composition, coverage, charge or temperature; return an out-of-domain flag to the caller.

Downstream sampling must reproduce MD moments and correlations. Round-trip test by drawing a large synthetic ensemble from the exported kernel and comparing probability, yield, energy, angle, products, momentum/energy balance, and tails to held-out raw trajectories.

Distinguish reset-surface and cumulative-dose experiments. Reset-surface trajectories isolate a conditional single-impact response at fixed state $\chi$. They are ideal for kernels but do not predict how the surface evolves. Cumulative bombardment updates composition, coverage, roughness, implantation and damage after each impact and can reveal steady state; its artificial numerical flux and finite reservoir demand care.

If the simulation area is $A$ and $N_{imp}$ impacts have occurred, represented fluence is

Physical time is $t=\Phi/\Gamma$ for incident flux $\Gamma$, but MD usually applies impacts many orders faster than experiments. The lattice may be equilibrated between shots, yet slow diffusion, desorption, adsorption and radical arrival are not automatically restored. Calling the sequence a real flux experiment without time-scale bridging is incorrect.

Use controlled protocols: reset to a sampled fixed state for conditional kernels; alternate impact MD with kMC for slow surface evolution; insert statistically justified relaxation/reservoir steps; or index kernels by measured/modelled state and let the feature model evolve it. Compare ordering and relaxation assumptions.

Cumulative slabs also deplete or enrich because products leave while incoming feed is limited. Maintain elemental inventory, surface recession reference, density, charge and temperature. Add/remove bottom layers or remap the slab only with conservative procedures. A 5 nm slab cannot supply an indefinitely etched bulk.

MD and surface kMC divide fast dynamics from rare time. MD resolves ballistic collision cascades, prompt reactions, product recoil, local thermalization and short correlated rearrangements. kMC advances activated adsorption/desorption, diffusion, reaction, restructuring and other events separated by long waits. Transition-state theory may supply a thermal rate

with barriers/prefactors from electronic structure or enhanced sampling, not ordinary short MD counts.

A coupled cycle may: sample state/site from kMC; extract an atomistic neighborhood; run an MD impact; classify products and map retained atoms/damage back; update kMC event catalog; then advance slow events until the next impact hazard. Mapping must conserve atoms, sites, coverage, damage and energy reservoirs. It also needs a policy for novel states that lack rates.

Avoid double counting prompt chemistry in both MD and kMC. Define a commitment time/surface separating fast impact outcome from slow events. Recrossings around that boundary require state definitions or milestoning. Validate coupled results against longer MD where accessible and limiting analytic cases.

Accelerated MD, metadynamics, hyperdynamics, parallel replica, transition path sampling and rare-event methods answer specialized questions, but their biased clocks/ensembles must be decoded correctly. Do not mix biased trajectories into an ordinary dynamical kernel without reweighting.

MD connects to feature Monte Carlo through boundary events. The feature model provides local material/state and incident $(s,E,\Omega)$; MD returns outcome/product samples or kernels, surface-state increments, deposited energy/momentum, implantation depth and uncertainty. The feature model transports reflected/products to their next hit and accumulates removal/deposition volume.

Material removal conversion needs an explicit counting unit and atomic/formula-unit volume. If $Y_m$ target units leave per incident particle and local incident flux is $\Gamma$, planar recession speed is

where $n_m$ is number density of the same unit. Mixed/porous/polymer states require state-dependent density and composition. Do not convert a silicon-atom yield using SiO$_2$ formula-unit density.

Return product ancestry so feature MC can distinguish direct reflection, sputter, etch product, and secondary reaction. Return outgoing particles only once at a declared crossing surface. Lattice heat and retained atoms belong to surface/material state, not new rays. Charged-current bookkeeping stays consistent with the model's neutralization assumption.

A library should span relevant material interfaces, mask/film states, temperature, energy, angle, species mixtures and dose variables. Sparse high-dimensional sampling needs designed experiments, adaptive acquisition, surrogate uncertainty, and physical constraints. A black-box neural surrogate without conservation or out-of-domain signaling is not a safe closure.

Finite-size and duration artifacts have recognizable signatures. Periodic-image correlations cause products/cascades to repeat with cell pitch. A thin slab shows bottom reflection or fixed-layer damage. Too-short runs classify slow products as adsorbed; too-long runs receive acoustic echoes. Too-strong damping lowers sputter and product energy. Too-small vacuum causes products to interact across periodic $z$ or return from a wall.

High instantaneous bombardment produces runaway heating, overlapping cascades, unnatural radical depletion, excessive damage and composition drift. Frequent global rethermalization erases real damage and retained heat. A single crystalline impact position produces channeling/site bias. An over-smoothed potential misses repulsion; a hard splice leaks energy.

Track diagnostic profiles versus depth/time: kinetic temperature by region, coordination, composition, defects, stress, deposited energy, maximum displacement, product crossing, and thermostat work. Plot them by event/dose and compare cell sizes. Visual inspection is useful for classifier debugging but cannot replace ledgers and convergence.

Surface location itself evolves. Define it using density/composition/height field or atomistic surface detection, then use the same definition for incidence angle, implantation depth, and escape plane. A fixed original plane gives misleading depths after recession/deposition.

Rough surfaces require local normal at the interaction scale. Atomic normals fluctuate; mesoscopic feature normals omit atomic corrugation. State whether angle-conditioned kernels average over equilibrium corrugation about a macroscopic normal. Avoid feeding atom-scale normal noise into a continuum ray model twice.

Verification asks whether the implemented model solves its stated equations. Test the integrator on harmonic/Kepler or known lattice systems; demonstrate NVE energy-drift order with timestep; verify force/energy consistency by finite differences; test neighbor-list skin/update under high-speed collisions; check periodic and nonperiodic boundaries; and reproduce temperature/velocity distributions after equilibration.

For the potential, reproduce published/fitted properties with the exact runtime settings and test held-out structures, forces, reaction paths, repulsive collisions, molecules and surfaces. Confirm unit conversions and element ordering. Compare energy/forces across code implementations for a small frozen configuration where possible.

For event analysis, build synthetic trajectories with known adsorption, crossing, return, fragmentation, mixed products, periodic wrap and implantation. Every atom must be counted once. Test classifier sensitivity to bond/distance/persistence and escape-plane position. Verify weighted and unweighted estimators on controlled random samples.

For the kernel, integrate normalization/multiplicity, sample it back, and close atoms, momentum and energy within declared reservoirs. Confirm interpolation at grid nodes, within cells, and at the domain boundary. Deliberately query out-of-domain inputs and require an explicit refusal/fallback.

Validation proceeds from elementary properties to process observables. First validate potential chemistry and surface structure. Then compare single-species beam experiments where energy, angle, material, temperature, coverage and dose are controlled: reflection, sputter/etch yield, threshold, product identity, angular/energy distributions, implantation/range and damage. Use multiple observables so compensating errors cannot hide.

Plasma exposure adds mixed radicals/ions, charging, photons, changing coverage and thermal/transport effects. Reproduce incident distributions with a sheath/diagnostic model and compare dose-dependent composition, film thickness, mass loss, product spectroscopy, roughness and etch yield. A match to one steady yield after tuning does not validate the product kernel.

Split uncertainty into potential/model form, reference electronic structure, initial surface ensemble, incident distribution, finite cell/time, classifier, sampling, experimental measurement, and multiscale mapping. Propagate influential components to final profile metrics. Do not hide model-form uncertainty inside narrow trajectory confidence intervals.

Calibration may update a small set of physically interpretable parameters or discrepancy terms using training experiments. Preserve held-out materials/energies/angles/recipes for validation. Avoid tuning the potential and surface kernel independently to the same final contour; it destroys identifiability and provenance.

Performance is an ensemble problem. Impact trajectories are naturally parallel across surface replicas and conditions. Use GPU/domain decomposition only after verifying force/energy parity, neighbor handling, RNG independence and reduction precision. Reactive charge equilibration or long-range electrostatics may dominate runtime; profile before approximating.

Record completed qualified impacts per compute-hour, including failed/out-of-domain trajectories. Checkpoint cumulative-dose runs with atomic state, velocities, integrator/thermostat state, RNG, inventories and event history.

Reproducibility requires input decks, initial configurations/ensembles, potential files and hashes/licenses, electronic-reference settings, software/compiler/device versions, seeds, unit system, timestep/protocol, raw trajectories or lossless event records, classifier code/version, convergence notebooks, and kernel schema. A paper citation to a force-field family is insufficient.