Charge Collection After Heavy-Ion Strikes: A Material Parameter Sweep in a 2001 TCAD Solver

Charge Collection After Heavy-Ion Strikes

A material parameter sweep in a 2001 TCAD solver

This is a follow-up to an earlier note about bringing ISE TCAD 7.0 — a
semiconductor device simulator from 2001 — back to life on an ARM Mac under
emulation. That note covered the infrastructure. This one covers what came out
of actually using it.

Short version: the simulation results are modest, and the most useful thing the
study produced was a lesson about how easy it is to fool yourself with transient
integration settings.

The setup

A two-dimensional silicon diode with an abrupt p-n junction. The mesh is
deliberately small — 162 vertices, 282 triangles — so that a single run
costs minutes rather than hours.

device      3.0 × 1.0 µm
anode       1.0 µm wide, top face (partial contact)
cathode     full bottom face
junction    at mid-depth
bias        −2 V reverse

The anode covers only a third of the top face. This turns out to matter.

Ion strikes are modelled with the solver's built-in heavy-ion source:

HeavyIon (
    Direction = (0,1)
    Location  = (0.0,-0.5)
    Time      = 100e-12
    Length    = 1e-4        ! cm, i.e. 1 µm
    Wt_hi     = 0.05e-4     ! cm, i.e. 0.05 µm track radius
    LET_f     = 0.2         ! pC/µm
    Gaussian
    PicoCoulomb
)

Two units traps worth flagging, because both cost me a failed run: Length and
Wt_hi are centimetres, not microns, and LET_f defaults to pairs/cm³
unless you add the PicoCoulomb keyword. Get either wrong and the model
silently generates almost nothing — no warning, no error, just a flat current
trace.

Validation 1: does the ion source work?

Baseline run on plain silicon. The track deposits 0.2 pC/µm over a 1 µm path,
so 0.2 pC total.

Anode current vs time, log scale

The current jumps twelve orders of magnitude and decays within a couple of
hundred picoseconds. Eighty-eight per cent of the deposited charge reaches the
anode; the rest recombines. Both the shape and the magnitude are what you would
expect, which is the point of running a baseline.

Validation 2: can we invent materials?

The solver loads material parameters from a plain-text library, and a
simulation can override any of them through a small file:

Material = "Silicon" {
    Bandgap {
        Eg0 = 1.40000   # [eV]
        ...
    }
}

Only the blocks you list are replaced; everything else falls back to the
library. That makes hypothetical materials easy: take silicon, change the
numbers you care about, leave the rest alone.

To check that this actually changes the physics — and changes it by the right
amount — the same diode was run twice with nothing different but the bandgap.

I-V curves for two bandgaps, log scale

At low bias the current is proportional to nᵢ², which scales as exp(−Eg/kT).
For ΔEg = 0.2304 eV at 300 K that predicts a ratio of

exp(0.2304 / 0.02585) = 7420

against a measured 7260 — agreement within 2 %. At high bias the curves
converge, because there the current is limited by series resistance, which does
not care about the bandgap.

So the override mechanism works, and it works quantitatively.

One piece of physics you must not forget

The heavy-ion source specifies deposited charge per micron, not deposited
energy. Sweep the bandgap while holding LET_f fixed and you are implicitly
claiming that a wide-gap material generates as many electron-hole pairs per unit
path as a narrow-gap one. It does not.

The energy needed per pair follows Klein's rule:

E_pair ≈ 2.8 · Eg + 0.6   [eV]

so at fixed deposited energy the pair count scales as 1/(2.8·Eg + 0.6). Across
the range swept here that is a factor of 2.67 fewer carriers at Eg = 3.3 eV
than in silicon. Every run below rescales LET_f accordingly.

Skip this and wide-bandgap materials look artificially good.

The sweep

Nine points: bandgap × electron mobility, everything else held fixed —
geometry, mesh, bias, strike position, deposited energy.

Eg   ∈ {1.1, 2.2, 3.3}  eV
µn   ∈ {200, 900, 3000} cm²/(V·s)
µp   = µn / 3

Each point takes about six minutes. The full sweep ran unattended in under an
hour.

Collection efficiency vs mobility, one curve per bandgap

Collection efficiency, %

Collected charge, pC

What the two parameters do

They act on different observables, and that separation is the main result:

• Bandgap sets the absolute signal. Collected charge falls 2.73× going from
  1.1 to 3.3 eV, almost entirely through Klein's rule. Efficiency barely moves.
• Mobility sets the efficiency. Five to ten percentage points across a
  fifteen-fold mobility range, at every bandgap tested. Absolute charge barely
  moves.

If you are fitting parameters to match a measured pulse height, the bandgap is
what you need to know accurately. If you are matching collected fraction, it
is the mobility. Effort spent on the wrong one is wasted.

The counter-intuitive part

Collection efficiency decreases as mobility rises. The naive expectation is
the opposite: faster carriers should escape before recombining.

I do not have an established mechanism. The leading candidate is diffusive
spreading. The diffusion coefficient follows mobility through the Einstein
relation, D = µkT/q, so fifteen times the mobility means fifteen times the
diffusivity. The track plasma spreads sideways as √(4Dt), and with a 1 µm anode
on a 3 µm wide device, a larger share of the carriers wanders out of the
collecting region and recombines.

That hypothesis is cheap to test: widen the anode to the full face and the
effect should shrink. Not done here.

The part that actually cost a day

One grid point — Eg = 3.3 eV, µ = 200, the corner where the plasma lives longest
and is extracted slowest — returned four different answers under four
integration schemes:

The correct value is 95.9 %, which makes that corner the best point on the
map. The first attempt had it as the worst by a wide margin.

Two things follow.

The integration scheme can outweigh the thing you are measuring. The spread
across those four artefacts is larger than the spread between many published
parameter sets for real materials. A material comparison run on a single fixed
time scheme can be dominated by the scheme rather than by the materials.

The cost of a correct answer is not uniform across parameter space. Eight of
nine points converged on the first scheme. The ninth needed four attempts — and
it was not a random point but the physically stiffest corner. Sweeps run at
uniform settings are therefore systematically less accurate exactly where the
physics is most interesting.

Neither observation is deep. Both are easy to fall for: the first artefact
produced a saddle-shaped map with a trend reversal, which looked like a finding.
It was a truncated tail.

The check that catches it

Cheap and sufficient: for every run, report the residual current at the end of
the window as a fraction of the peak.

        µ=200      µ=900     µ=3000
Eg=1.1  1.4e-05    1.2e-12   4.0e-13
Eg=2.2  5.1e-07    2.2e-14   3.0e-15
Eg=3.3  3.3e-01 ←  9.6e-15   1.9e-16

The one bad number stands out immediately. Two points were additionally
re-run with a sixty-times-longer window and reproduced their original values to
a tenth of a per cent, which is what makes the rest of the table trustworthy.

Limits of this study

Stated plainly, because they are not small:

• The grid is 3×3. Enough for trends, not for a response surface.
• Parameters were swept independently. In real materials bandgap,
  permittivity and mobility are correlated; some of these combinations cannot
  exist.
• Ranges were chosen by hand, not taken from literature spreads.
• Hole mobility was set to one third of the electron mobility with no
  justification beyond plausibility.
• One geometry. The mobility effect is very likely tied to the partial
  anode, so nothing here transfers to another device unchecked.
• No experimental validation, and none possible with this setup.

This is a pilot, not a result to build on.

Reproducibility

Everything is scripted from a Jupyter notebook on the host machine; the solver
runs in the virtual machine and is driven over a network shell. Total wall time
for the nine-point sweep was under an hour of unattended computation.

The solver itself is fast: two Newton iterations per timestep, 0.13 s each, on
this mesh. Long runs are cheap once they are prevented from dying — which, as
the previous note describes, was its own separate problem.

Tooling written along the way:

• a DF-ISE parser, because the solver's output format encodes elements by edges
  rather than vertices and splits fields across regions
• a plotting layer reproducing the original visualiser's conventions
• a parameter-file generator for material sweeps
• a network shell client with prompt synchronisation

Would I take this further?

Not as research. The findings are thin, the mechanism is unexplained, and the
methodological observations — check your integration window, check your
timestep — are lessons rather than discoveries.

As an exercise it was worth it. The stand works, the pipeline is reusable, and a
solver old enough to vote turned out to be perfectly capable of producing
quantitatively correct physics, provided you are careful about units, about time
integration, and about not believing the first plot that looks interesting.