Validation & Verification (V&V)

This simulator uses a reduced fluid formulation for fast relative-trend analysis of RF/CCP plasma. It estimates qualitative and semi-quantitative trends in plasma density, ion flux, RF impedance, sheath response, and parameter sensitivity. It is not a replacement for fully validated PIC/MCC, commercial multi-fluid simulation, or experiment.

• Recipe trend comparison (A vs B under the same gas conditions)
• Pressure / power / gap / frequency sensitivity analysis
• Relative plasma-density trend estimation
• Relative ion-flux trend estimation
• RF impedance and matching tendency exploration
• Educational visualization and physical intuition building
• Early-stage process intuition (design-space screening)

Do not rely on this simulator for the following purposes — results may be misleading.

• Absolute density values without calibration
• Final process qualification or production sign-off
• Detailed radical chemistry analysis
• Low-pressure nonlocal kinetic analysis
• Stochastic sheath-heating analysis
• High-frequency electromagnetic wave effects
• Hardware design sign-off
• Replacement for experiment, PIC-MCC, or validated commercial solvers

Uncalibrated electron density uncertainty is approximately 2–5× of the reference. Absolute ion flux is accurate to roughly one order of magnitude. Calibration with experimental data improves relative trend accuracy near the calibration operating point, but absolute accuracy remains system-dependent. For applications requiring absolute confidence over a wide parameter range, use full PIC/MCC or a validated multi-fluid commercial solver.

The upgrades below fall into two groups. (a) Always-on — applied to every run with no opt-in. (b) Flag-gated — controlled via environment variables; ION_SOLVER, RF_SHEATH, and TE_PDE default to ON ("0" to disable), while only GUMMEL defaults to OFF (bit-exact when OFF). This remains a trend-level reduced fluid model and does not guarantee absolute accuracy.

(a) Always-on (enabled by default)

(b) Flag-gated (ION_SOLVER, RF_SHEATH, TE_PDE default ON; only GUMMEL default OFF)

ION_SOLVER, RF_SHEATH, and TE_PDE are now ON by default and can be disabled by setting the environment variable to "0". Only the GUMMEL outer loop remains OFF by default (opt-in); when OFF, results are bit-exact with the existing pipeline.

1. Run a baseline Ar simulation at the documented GEC-cell conditions (Hargis et al. 1994).
2. Compare computed electron density and ion-flux trends against published PIC/MCC or fluid benchmark results at the same conditions.
3. Confirm that trend directions (increasing/decreasing) agree across voltage, pressure, and gap sweeps.
4. If available, derive a calibration factor from experimental data (Langmuir probe, VI probe, matching network).
5. After calibration, trust absolute numbers only near the closest operating point.
6. Always review the Scope & Limitations section before drawing conclusions.

• GEC reference cell (Hargis et al. 1994) — the well-documented parallel-plate 13.56 MHz Ar CCP benchmark experiment.
• Ar CCP electron density and ion flux — compare directly against literature values (density ~10⁹–10¹¹ cm⁻³, ion flux ~10¹⁵–10¹⁶ m⁻²s⁻¹).
• Pressure / power / gap / frequency sweep trends — check directional agreement against Boeuf & Pitchford (1995) and Surendra & Dalvie (1993).
• VI-probe and matching-network impedance — compare RF impedance trends against lab measurements.
• COMSOL Plasma Module — a validated commercial multi-fluid solver for trend comparison.
• Zapdos/MOOSE — open-source MOOSE-based plasma finite-element code for independent cross-check.

The table anchors each model output to published-literature values. Each row lists a preset deep-link (loads directly into the simulator), the model output to read, the literature anchor, and the pass criterion. These are order-of-magnitude and directional checks, not absolute validation.

Explicit exclusions (out of scope for literature validation)

The following are deliberately out of scope for literature validation (mirrors physics-model.md §12).

• Process outcomes of every kind — etch rate, deposition rate, selectivity, profile (radical chemistry and surface reactions are not modeled; plasma quantities only, and no relative proxy either)
• EEDF / non-Maxwellian detail
• Space-charge sheath potential profiles (deliberately rejected — unresolved λ_De)
• Absolute IEDF eV scale (power absorption / |E| are now in physical units W/m³, V/m, power-balance normalized to the requested power)
• IEDF bimodal structure (no RF-cycle resolution)
• Stochastic sheath heating (ohmic Drude only)
• Absolute impedance (ohms)
• Full CFD neutral flow (now a showerhead→pump potential flow — face-selective emission + jet-column directionality — but not Navier–Stokes CFD)

The matrix lists the expected trend and pass/fail criterion for each representative benchmark case. These are directional and rough-magnitude (ratio / order-of-magnitude) checks, not absolute-value accuracy. It complements B13 (literature anchors V1–V10). Status chips: verified (green) = the full criterion is asserted by an automated test in tests/test_validation_matrix.py; partial (blue) = the direction is CI-tested but one clause of the criterion is not; expected (neutral) = literature-grounded but not yet a regression check; todo (amber) = no coverage.

The explicit out-of-scope list (will not match) is kept in one place, → What we don't claim above.

References below are grouped by model component. References already cited in the sections above are not duplicated.

Drift-diffusion / fluid model (general)

• Lieberman & Lichtenberg, Principles of Plasma Discharges and Materials Processing, 2nd ed., Wiley 2005 — drift-diffusion, global model, sheath (already cited in Section 1)
• Boeuf & Pitchford, Phys. Rev. E 51, 1376 (1995) — RF fluid CCP model, Poisson-coupled fluid + power absorption (already cited in Section 1)
• Scharfetter & Gummel, IEEE Trans. Electron Dev. 16, 64 (1969) — exponential-fit flux discretization (already cited in Section 3)

Electron rate-coefficient LUT / EEDF (A1)

• Hagelaar & Pitchford, Plasma Sources Sci. Technol. 14, 722 (2005) — BOLSIG+ Boltzmann solver for electron transport coefficients and rate coefficients
• BOLSIG+ official site (LAPLACE, Toulouse) — free electron Boltzmann solver for generating LUT data
• LXCat cross-section database — source data for A1 LUT / BOLSIG+ inputs (electron-neutral cross sections)

0D electron energy / particle balance (A2)

• Lee & Lieberman, J. Vac. Sci. Technol. A 13, 368 (1995) — global model 0D energy and particle balance for high-density plasmas

Sheath / Bohm criterion / RF sheath (A3)

• Riemann, J. Phys. D 24, 493 (1991) — Bohm criterion and sheath formation (already cited in Section 4)
• Lieberman, IEEE Trans. Plasma Sci. 16, 638 (1988) — analytical RF sheath model (already cited in Section 4) [verify IEEE Xplore]
• Godyak & Sternberg, Phys. Rev. A 42, 2299 (1990) — dynamic model of electrode sheaths in symmetrically driven RF discharges

GEC benchmark cell (B8)

• Hargis et al., Rev. Sci. Instrum. 65, 140 (1994) — GEC RF reference cell: defined parallel-plate 13.56 MHz system for benchmark studies

Comparison / open-source solvers (B8)

• COMSOL Plasma Module — validated commercial multi-fluid plasma solver for benchmark comparison
• Zapdos (shannon-lab / MOOSE) — free and open-source MOOSE-based finite-element plasma simulation code

When the magnetic field (B) input is active, the simulator applies a trend-level magnetized-electron confinement model. Core physics: the static B field magnetizes electrons, forming the Hall parameter β = ω_ce/ν_m (ω_ce = eB/m_e: electron cyclotron frequency, ν_m: electron-neutral collision frequency). Large β reduces cross-field electron transport by the factor 1/(1+β²), applied directionally as aniso_r and aniso_z. The result is enhanced electron confinement and a density increase trend captured by mag_gain.

This magnetic field model does NOT include: full Hall-tensor anisotropic transport / ExB-drift rotation / electron cyclotron resonance (ECR) heating / magnetized-sheath physics. It is a trend-oriented approximation based on directional Hall-parameter anisotropic transport + a lightweight ExB net-trend gain. Do NOT use it as an absolute magnetized-plasma result.

• Lieberman & Lichtenberg 2005, Ch. 4 — electron transport in magnetic fields, Hall parameter
• Raizer 1991, Ch. 9 — magnetized electron transport and drift-diffusion in magnetic fields

The table below summarises the relative trend directions this reduced-fluid model reproduces, verified by parameter sweeps. Arrows (↑ ↓ ~) indicate which way outputs move — they are not calibrated absolute magnitudes. For absolute use, calibration against experimental data is required.

↑ = increases, ↓ = decreases, ~ = approximately flat. Magnitudes are trend-level; calibration against measurements is needed for absolute use.

Te(r,z) map (fields.te_eV)

The fields.te_eV heatmap is, by default, a heuristic energy map mean-matched to the calibrated scalar Te (a geometric blend of the heuristic and the 0D particle balance, clamped to [1,9] eV). With PLASMACCP_TE_PDE, the spatial shape instead comes from a screened-diffusion PDE heated by the RF absorption map.

nₑ observable (fields.ne)

The displayed fields.ne = v/(v+n_sat) ∈ [0,1] is a saturating transfer (a per-request fixed curve so A/B legends stay comparable). The absolute density scale lives in density_ref / the 0D block.

Electrode gap

Gap-derived quantities (the Te gap factor, the Child-Langmuir clamp, and Λ=gap/π) use the grid-derived discharge gap — the median plasma-cell span between the powered and ground blocks. z_max is used only as a fallback.

Calibration knobs (calibration block)

Neutral defaults (1.0 / 0.0) = bit-exact with prior results.

B1–B14 above describe what is modeled; this section documents how it is actually computed — the grid, the discretization, the linear-algebra solvers underneath, and the verification gates layered on top.

Numerical honesty: the goal is trend direction, not absolute value. CI pins on every commit — parity of every vectorized path against its scalar reference, off=bit-exact regression for every new flag, directional assertions (power↑→density↑, etc.), an r-alternation (checkerboard) guard, and the M1–M8 benchmark suite (must pass unmodified).

Bit-exact escape hatches: PLASMACCP_PLANAR_RADIAL_OPERATOR=1 (reproduces the old planar operator), PLASMACCP_TG_FEEDBACK=0, calibration neutral values (1.0/0.0), and RF_SHEATH="0" (Child-law fallback only). The native C assembler stays disabled by default in prod — it carries a checkerboard M-matrix bug at high cell aspect ratio and, since BC2, is also physically stale (still planar); a lone ENABLE=1 is ignored (with a warning) by a consistency guard unless paired with PLANAR_RADIAL_OPERATOR=1.

The V1 engine has two ways to set the absolute electron density. Standard mode (today's default) uses the calibrated empirical density scale as-is. Physical mode (beta) instead derives the density level directly from an inter-electrode wall power balance — absorbed power divided by the energy ionization loses at the walls. Both run on the same 2-D density-distribution solver; physical mode only switches how its absolute scale is set.

Mechanism

Density is solved from a per-face power balance n = P_abs / (e·u_B·(E_c(Te)+2Te+V_sheath)·A_eff). The ionization source sits at the driven sheaths — the powered electrode and the facing electrode (including the conductor behind a wafer sitting on it) — smoothed spatially over the electron energy relaxation length λ_ε (Tsendin 2010). The sheath voltage is re-solved self-consistently (iterated to convergence) against that same power balance. Finally a Bohm self-consistency check (R1) compares the ion current density expected from the solved bulk density and Te against the actually reported ion current density, per channel (powered electrode, wall) — outside a factor-of-κ=4 band it blends back toward the empirical scale in log space, proportionally to how far outside the band it fell.

R1 checks only self-consistency (whether a channel's own density and Te agree with that channel's own current) — it is not checked against an external literature value, so passing R1 is not evidence that the absolute density is correct.

Te is also solved from a particle balance for every gas and mixture, using the three-term Lee–Lieberman h-factor and a per-gas ion mobility (measured for Ar, N2, O2, H2 and He; estimated for the rest, and the run notes say which). The Te shown is the one the density calculation actually used. Multi-frequency requests get the physical density too: the tones drive one plasma on their combined power, and the power balance's sheath voltage is solved per tone at its own power and summed. The attachment loss applies only to gases that actually attach.

Gas coefficients (ionization rate k_iz, energy cost per ionization E_c) rest on different evidence per gas: Ar, He, N2, O2, H2, SiH4, CF4, Cl2 and N2O have a dedicated cited fit (N2O's is integrated directly from Song et al. 2019 cross-sections); NH3, SF6, CH4 and C2H2 use an integral of NIST's public BEB ionization cross-sections; Ne, Kr, Xe, NF3, C4F8, CHF3 and HBr reuse a chemically similar gas's fit shape, re-anchored only by that gas's own cited ionization potential (not a newly fitted coefficient). A gas with no entry at all still returns a result — it falls back to argon's coefficients, and the run notes say so explicitly. A mixture's energy cost blends each component's own value weighted by how much it actually ionizes, not by its mole fraction alone.

For multi-frequency requests, the process monitor shows each tone's sheath voltage and ion power separately (Lieberman's dual-frequency ion-energy formula applies when there are exactly two tones).

The three badge states

Validation status — honestly

• Completion (wide range): a 950-cell sweep (475 standard, 475 physical) across 26 gases/mixtures (including 11 newly added) combining randomized (LHC) and one-factor-at-a-time conditions completed in full — zero exceptions or refusals, every reported number finite, runtime within 3x standard (median 0.78x). 85.9% kept the physical density outright (wall power balance in the table above); 14.1% fell all the way back to the empirical anchor because the self-consistency check (R1) landed outside its band — this counts as degraded, not failed, and mostly occurred at extreme randomized combinations. R1 is only a self-consistency check, so this does not mean the absolute values are right. (The earlier, narrower 268-cell sweep — 10 gases x pressure only — remains the basis for the absolute-value, shape, and power/pressure-dependence items below.)
• Absolute value — not yet validated: we compared against the Turner et al. (2013) He PIC benchmark, a public reference with absorbed power explicitly defined. In a flat-plate chamber with the gap matched to Turner's 67 mm, the centre density came out at 0.85–4.8x Turner's value in physical mode (0.85–1.5x at three of the four pressures) and 0.06–0.62x in standard mode. In a small research cell (25 mm gap, earlier measurement), physical mode gave 0.11–1.32x and standard mode 0.11–0.62x. The true value appears to lie between the two modes — neither is a validated absolute number.
• Profile shape: in an earlier round, physical mode's density profile had a non-physical dip at mid-gap centre (from a duplicated volume-loss term). After removing it, all four of Turner's benchmark cases come out centre-maximum with a flat top, matching the shape of Turner Fig. 7–10.
• Power dependence: at GEC conditions (26.7–53.3 Pa), physical mode's power exponent is 0.85–0.90, close to the measured value (Land et al. 2009, ≈0.98–1.08). Standard mode's exponent is 0.086 — essentially flat over that power range.
• Pressure dependence: against the Turner benchmark (30→1000 mTorr, ×18.4), physical mode gives ×6.6–14.5 and standard mode ×2.6–3.3 — physical is closer, both are too shallow. At GEC conditions (26.7→53.3 Pa) both modes rise with pressure (physical exponent about +1.1). Land et al. (2009) measured a fall, but their power is apparent power (V·I without phase) and the authors themselves note the result differs from the literature. Public references with (nearly) fixed absorbed power — Turner 2013 PIC, Schulenberg et al. 2021, the NIST GEC special issue (Hopkins, Overzet 1995) — all show a rise. The heating-mode transition (Godyak & Piejak 1990, about 0.4–0.5 Torr) is absent in both modes: V1 does not model it.

When to use it, and when not to

• Use it for: exploring how density responds when power is varied — physical mode's power exponent (~0.9) tracks measurement far better than standard mode's (~0.09).
• Do not use it for: anything that needs an absolute density number (process qualification, tool-spec comparisons) — against a benchmark with defined absorbed power, physical mode is up to about 5x too high depending on conditions, and its ion energy (sheath voltage) is far too low. The exact pressure exponent and the heating-mode transition region are unvalidated.
• Standard mode remains the default; physical mode is an opt-in beta. Always read it together with the density-source badge (table above) — when it shows "blended" or "empirical anchor," that result's absolute density came from the same path as standard mode, not from the physical calculation.
• Known limits — stated explicitly: (1) Electronegative gases (Cl2, CF4, SF6, etc.) are still solved as electropositive — attachment feeds into the rate coefficients, but negative ions are not tracked as a separate charge carrier. (2) In the 950-cell sweep above, physical density falls all the way back to the empirical anchor in 14.1% of cases — in those conditions physical mode produces the same numbers as standard mode. (3) He and NH3 have no openly licensable measured cross-section data we could obtain, so terms beyond the ionization rate (He's excitation/elastic terms, NH3's energy cost and attachment) are still proxies or borrowed from another gas.

• Turner, Hutchinson, Kushner et al., Phys. Plasmas 20, 013507 (2013) — arXiv:1211.5246 (open access); multi-code PIC-MCC benchmark, He, absorbed power explicitly defined
• Lieberman & Lichtenberg, Principles of Plasma Discharges and Materials Processing, 2nd ed., Wiley 2005 — wall power balance / global model (already cited above)
• Tsendin, Phys.-Usp. 53, 133 (2010) — electron energy relaxation length λ_ε, used to smooth the fixed ionization source

The Process summary tab reports wafer-level trend indicators. The subsections below explain what each metric means, why a process engineer cares, and how the simulator computes it.

Radial ion-flux nonuniformity

Measures the radial flatness of the ion-flux profile. A larger value means the ion flux arriving at the wafer centre and at its edge differ by that much. The simulator computes nonuniformity (%) = (max − min) / (2 × mean) × 100 over the radial profile; lower is better. If the offered wafer interval contains a zero or non-finite sample, the scalar is hidden and used/total column coverage is shown because missing surface support cannot be distinguished from a physical zero. Real-tool uniformity also depends on gas-flow patterns, chuck temperature, and focus-ring condition — factors not fully captured here.

• Lieberman & Lichtenberg, Principles of Plasma Discharges and Materials Processing, 2nd ed., Wiley 2005 — Ch. 10–11 (uniformity, transport)
• Boeuf & Pitchford, Phys. Rev. E 51, 1376 (1995) — 2D CCP fluid model with radial uniformity analysis

Ion energy / IEDF

The energy of ions impinging on the substrate is set by the sheath voltage drop, and it fixes how much kinetic energy each ion delivers to the surface. The shape and mean of the ion energy distribution function (IEDF) describe that distribution. The simulator computes a time-averaged ion energy proxy based on Child-Langmuir scaling at each radial position — this is a relative trend indicator, not an absolute eV value. In dual-frequency configurations, the high frequency (HF) primarily controls ion flux while the low frequency (LF) primarily controls ion energy (Boyle et al. 2004). Real IEDFs show a characteristic double-peak structure arising from RF sheath phase-averaging (Godyak et al. 1991); this simulator does not resolve that structure.

• Lieberman & Lichtenberg 2005, Ch. 11 — RF sheath and IEDF theory
• Boyle, Ellingboe & Turner, J. Phys. D 37, 697 (2004) — HF controls ion flux, LF controls ion energy in dual-frequency CCP
• Godyak, Piejak & Alexandrovich, J. Appl. Phys. 69, 3455 (1991) — measured ion energy distributions and flux in RF discharges
• Lieberman, IEEE Trans. Plasma Sci. 16, 638 (1988) — analytical RF sheath model [verify IEEE Xplore] [verify]

Sheath thickness

The sheath is the positive space-charge region near electrodes where ions are accelerated toward the substrate. Sheath thickness is governed by the applied voltage, plasma density, and frequency (Child-Langmuir law), and directly influences ion bombardment energy and angular distribution. A thicker sheath generally corresponds to a longer ion acceleration path and higher ion energy. The simulator estimates sheath thickness at each radial position using the Bohm velocity criterion and Child-Langmuir scaling — interpret as a trend signal, not an absolute measurement.

• Lieberman, IEEE Trans. Plasma Sci. 16, 638 (1988) — analytical RF sheath model [verify IEEE Xplore] [verify]
• Riemann, J. Phys. D 24, 493 (1991) — Bohm criterion and sheath formation
• Lieberman & Lichtenberg 2005, Ch. 2 & 11 — sheath physics and Child-Langmuir scaling

Ion flux–energy index

A single dimensionless index built from two plasma quantities: ion flux and ion energy. At each radial position the simulator forms I(r) = ion flux × an ion-energy weight — the threshold sputter yield √(max(E − E_th, 0)) when an absolute energy scale exists, otherwise √(relative energy) or no weight at all — and normalizes to an area-weighted mean of 1. The index describes the radial distribution of ion bombardment arriving at the surface.

The index is not a process outcome. This simulator does not predict etch rate, deposition rate, selectivity, or profile — those are dominated by radical chemistry and surface reactions this solver does not model. The physical sputtering that the index's energy weight describes is only one term of any real process, and dividing the index by anything does not recover a process outcome.

• Coburn & Winters, J. Appl. Phys. 50, 3189 (1979) — ion/neutral synergy at a surface: the measurement showing why an ion-only index cannot stand in for a surface outcome [verify DOI] [verify]
• Lieberman & Lichtenberg 2005, Ch. 11 — RF sheath physics setting the ion flux and ion energy at the electrode