Depletion handling: clamp bistability and the hidden 0 V floor

> [!CAUTION]
> **Correction (2026-07-21):** the central ground-truth claim of this issue — that the full-impurity *linear* solution at 1.05 × V_d shows no interior extremum, making the `6d576c13` "undepleted" result a spurious fixed point — is **wrong and retracted**. Direct linear solves (full impurities, `depletion_handling = false`) show the pocket above the point contact is **real ICPC pinch-off**, persisting to ≈ 1.08–1.12 × the estimated V_d; it is also reproduced in Cartesian coordinates. `6d576c13` reports genuine physics, and the current (ghost-anchored) code suppresses it — the ±0.5 % test brackets validate that artifact. The flawed original ground truth came from an adapted-weighting-potential superposition that is not faithful near the pocket. See the [correction comment](https://github.com/JuliaPhysics/SolidStateDetectors.jl/issues/618#issuecomment-5038538907) and discussion below for the corrected analysis and agreed direction (adopt `6d576c13`, re-anchor the depletion tests, initialize depletion handling from the converged linear solution).
>
> The parts of this issue that remain valid: the clamp's plateau fixed points and history dependence, the 0 V-floor fragility (including its gauge non-invariance), and the reformulation proposals.

Proposal generated as part of full AI review of SolidStateDetectors codebase:

From the 2026-07 review on the `improvements` branch; line numbers refer to that branch.
This documents the most consequential *algorithmic* finding of the review: `handle_depletion`
admits stable spurious "undepleted" states, and the thing currently preventing them in the
common cases is an accident — a never-written ghost value acting as a 0 V potential
floor at the cylindrical axis.

## How `handle_depletion` works

`src/PotentialCalculation/SuccessiveOverRelaxation/SOR.jl:57–83`, applied per grid
point inside the SOR update on CPU and GPU. Given the new potential proposal, the
six neighbor potentials, and the point's impurity term
`imp_contribution = q_eff_imp * volume_weight` (negative for p-type):

```julia
vmin = min(neighbor_potentials...); vmax = max(neighbor_potentials...)
neighbor_relevant_extremum = is_ptype ? vmin : vmax
imp_scale = (neighbor_relevant_extremum - new_potential) / imp_contribution
imp_scale = clamp(imp_scale, 0, 1)
if vmax == vmin imp_scale = 0 end
new_potential += imp_scale * imp_contribution
# overshoot clamp: never cross the neighbor extremum
if (is_ptype && new_potential < vmin) || (!is_ptype && new_potential > vmax)
new_potential = neighbor_relevant_extremum
end
```

The rationale (docstring, same file): in a depleted detector the potential is
monotone between contacts, so no interior point may become a local extremum. Where
the full impurity term would create one, the impurity density is scaled down
(`imp_scale < 1`) and the point is treated as undepleted; `is_depleted(point_types)`
reports `imp_scale == 1` everywhere.

## The accidental 0 V floor at r = 0

The cylindrical red-black layout is (z, φ, r). For an axis point the
radially-inward neighbor — slot 5 of the neighbor tuple — reads the sub-axis ghost
row, which **no boundary condition ever writes**: it stays at its initial value 0.
`r0_handling_depletion_handling`
(`src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl:236–250`)
keeps that slot in the tuple, so every axis point's extremum check runs against
`min(real neighbors…, 0)` — an implicit 0 V floor (ceiling, for n-type). The
Cartesian variant is a no-op: Cartesian grids have no such anchor anywhere.

This looks exactly like a stale-index bug, and the review first "fixed" it
(commit `6d576c13`: replace slot 5 by the outward neighbor). That fix made the
suite fail in a physically revealing way, which led to the analysis below and a
documented revert (`f87dc4ee`, plus a unit test pinning the tuple behavior in
`test/test_depletion.jl:7`).

## What the failed fix revealed: erosion fixed points

With the ghost masked, the inverted-coax example (`test_real_detectors.jl:26–37`,
estimated depletion voltage ≈ 1871 V) still reported **undepleted at 1.05× its
depletion voltage**. Ground truth says otherwise: because the no-depletion-handling
problem is linear, the true full-impurity solution at any bias is the superposition
`epot(V_b) = V_b · wpot + epot_ρ` (weighting potential of the bias contact, plus
the all-contacts-grounded impurity-only solution). Computing both on the refined
grid and scanning the superposition at 1.05× V_d shows no interior extremum — the
potential minimum sits on the 0 V contact. Depleted. The "undepleted" flag is a
spurious *fixed point of the iteration*, not a property of the physics:

* Wherever the local impurity term exceeds the curvature margin of the
neighborhood (|imp_contribution| ≳ vmax − vmin), the clamped update degenerates
to `new_potential = vmin` (p-type). The update becomes a **morphological
erosion**: each point copies the minimum of its neighbors.
* Erosion's stable states are flat plateaus pinned at the lowest value on the
region's rim. Flatness then keeps `imp_scale ≈ 0` (and `vmax == vmin` forces
exactly 0), which keeps the impurity term out of the update, which keeps the
plateau flat: self-consistent, and it survives arbitrarily many iterations and
further grid refinement.
* The true solution (imp_scale = 1, curved potential) is *also* a fixed point.
Which basin the solver lands in depends on history — thin, nearly-flat channels
(the IVC axis channel above the point contact) get trapped in the spurious one.

The ghost 0 is what breaks this at the axis: `vmin` is 0 there regardless of any
plateau value, so axis points are never pinned to a plateau neighbor, keep
(essentially) their full impurity term, and follow the true Poisson update; correct
values then propagate radially outward and drain axis-adjacent plateaus.

A second attempted fix — forcing `imp_scale = 1` at the axis — cured 1.05× but
broke the opposite test: at 0.995× V_d the detector must report *undepleted*, and
the legitimate signature of that at the axis is precisely "the impurity term would
push the point below the 0 V floor". The floor plays both roles; the current
behavior (revert + in-place documentation) is the only variant of the three that
passes 0.995× → undepleted, 1.005× → depleted, 1.05× → depleted.

## Why the current state is still fragile

1. **The floor hard-codes 0 V.** It is exact only while some contact is grounded
at 0 V (true for every shipped example). A bias scheme with both contacts away
from 0 (split bias, shifted reference) makes 0 wrong in both directions:
legitimate potentials below 0 near a negative contact read as undepletion, and
the anchor sits at the wrong level for the extremum check.
2. **The anchor exists only on the cylindrical axis.** Cartesian simulations and
off-axis flat regions (thin gaps between segmented contacts, plateaus formed
mid-refinement) run the same clamp with no floor at all; the erosion fixed
points are reachable there. (No failing Cartesian example was constructed
during the review, but it is the identical code path minus the accident that
protects the cylindrical axis.)
3. **`vmax == vmin ⇒ imp_scale = 0`** declares exactly-flat neighborhoods
undepleted, locking constant-initialized regions until rounding noise breaks
the tie.

## Proposed reformulation

Ordered by effort; 1 is a contained change, 2–3 are the real fix.

1. **Replace the implicit 0 by the applied-potential extrema.** For a depleted
detector the maximum principle bounds the potential by the applied contact
potentials, so `floor = minimum(applied potentials)`,
`ceiling = maximum(applied potentials)` are the correct generalization of the
ghost's role — and they can be applied uniformly at *every* point (extend the
extremum check by the bound), not just the axis. This fixes fragility 1,
removes the cylindrical/Cartesian asymmetry (2), and anchors every channel
whose erosion would drag it to a bound. It does **not** remove the bistability
for plateaus strictly inside the bounds.
2. **Make the spurious plateaus non-stationary.** Options that keep the current
structure: (a) when the overshoot clamp fires, mix a small fraction of the
*unclamped* full-impurity proposal into the update, so exact plateaus drift
instead of self-reproducing; (b) periodically re-anneal — reset
`imp_scale = 1` in flagged regions every N iterations and let the clamp
re-derive it; genuine undepleted regions re-form, spurious ones converge away;
(c) after convergence, re-evaluate flagged regions once with full impurities
and check for an actual interior extremum (the superposition criterion above),
warning or re-iterating on mismatch.
3. **Longer term:** what is being solved is a free-boundary/complementarity
problem (either full impurity and the potential respects the bounds, or the
potential sits on the bound and the charge is reduced). That is the classical
obstacle problem, and **projected SOR** (Cryer) solves its LCP form with a
convergence guarantee — a natural fit since the solver is already SOR. This
would replace the heuristic clamp entirely and make the depleted solution the
unique attractor.

Any change here alters physics-facing behavior, so it should be validated the way
the review validated the revert: the ±0.5 % depletion-voltage brackets for the
example detectors, plus superposition ground truth on refined grids, plus at least
one case with no 0 V contact and one Cartesian flat-channel case (which today have
no safety net).

## Reproducing the spurious state

*(Edited: the snippet originally posted here was flawed — without weighting
potentials at proper refinement, `estimate_depletion_voltage` returns ≈ 2156 V
instead of ≈ 1871 V and every variant reports depleted at all probe points.
Verified snippet below; see the comments for the full variant matrix.)*

With slot 5 masked in `r0_handling_depletion_handling` (the one-line change of
`6d576c13`), the settings of the "Inverted Coax" testset report **undepleted at
1.005× and even 1.05× V_d**; the current code reports depleted at both:

```julia
using SolidStateDetectors, Unitful
T = Float32
sim = Simulation{T}(SSD_examples[:InvertedCoax])
simulate!(sim, convergence_limit = 1e-6, refinement_limits = [0.2, 0.1, 0.05, 0.03, 0.02, 0.01], verbose = false)
deplV = estimate_depletion_voltage(sim, check_for_depletion = false) # 1871.4 V
id = SolidStateDetectors.determine_bias_voltage_contact_id(sim.detector)
for f in (0.995, 1.005, 1.05)
sim.detector = SolidStateDetector(sim.detector, contact_id = id, contact_potential = f * ustrip(deplV))
calculate_electric_potential!(sim, depletion_handling = true)
@show f, is_depleted(sim.point_types)
end
```

Relevant history on the branch: `6d576c13` (attempted fix), `f87dc4ee` (revert
with documented semantics), `test/test_depletion.jl:7` (tuple pin),
`SOR.jl:31–56` (docstring incl. the existing fixed-charge caveat, which this
analysis reinforces).

Start with src/PotentialCalculation/SuccessiveOverRelaxation/SOR.jl:31–83 and the CPU axis handling in CPU_innerloop.jl:236–250, then read test/test_depletion.jl:7 and the Inverted Coax case in test_real_detectors.jl:26–37. Validate any change against the ±0.5% depletion brackets, refined-grid superposition ground truth, a bias scheme without a 0 V contact, and a Cartesian flat-channel case; done means spurious plateau states are avoided without breaking legitimate undepleted results.

Written by the indexing model from the issue text.