Open-Source TCAD for the SkyWater Sky130 Process, Calibrated Against Measured Silicon from 300 K down to 77 K
Joonatan Alanampa — July 2026
(CI re-runs the stage 0–3 verification chain — solver-vs-theory and
solver-vs-DEVSIM physics checks — on every push, on a clean Linux machine.)
Repository = publication: every number below is produced by a script in this repo, every figure is committed, and every fitted quantity is explicitly separated from every predicted one.
We demonstrate a fully open-source technology-CAD (TCAD) workflow for the
SkyWater Sky130 130 nm process, running natively on a consumer Windows 11
machine, and calibrate it against three independent sets of measured
silicon: Google/SkyWater's published wafer characterization data at 300 K,
the same dataset's minimum-length (0.15 µm) device, and the UTA/Fermilab
liquid-nitrogen (77 K) dataset. A from-scratch 1D drift-diffusion solver is
first verified against closed-form theory and then against DEVSIM (agreement
to 0.1 % over five decades of current). A 2D nfet_01v8 model with a
field-dependent surface-mobility law reproduces the measured long-channel
transistor with a saturation-regime prediction error of ~6 %, and the
calibration transfers to the 0.15 µm short-channel device, reproducing the
full output-characteristic family to 1–9 % with series resistance modeled
and only a single on-current point fitted, and predicting DIBL within a
factor of 1.6. Because Sky130's doping profiles are proprietary, the
calibration itself measures them: the extracted channel doping is
N_A ≈ 3.9×10¹⁷ cm⁻³ (from the body effect), the short-channel device
requires an effective halo doping of 1.07×10¹⁸ cm⁻³ (from the −410 mV
threshold roll-off a halo-free model exhibits), and the source/drain
contact resistance extracts as R_c ≈ 380 Ω·µm. At 77 K, with dopant freeze-out (incomplete ionization) added to
Poisson's equation, the model matches the measured threshold voltage to 1 mV
and extracts a low-field mobility ratio µ₀(77 K)/µ₀(300 K) ≈ 6.7; its
predictions land in the same ~20–30 % error class as the dataset authors'
own isothermal BSIM4 compact models. Most notably, the ideal-physics
simulation quantifies the cryogenic subthreshold-swing saturation anomaly:
Boltzmann statistics demand 23 mV/dec at 77 K while the silicon measures
71 mV/dec — a 3× gap attributable to band-tail/interface-trap conduction
that standard drift-diffusion cannot capture. The calibrated TCAD
is exported as a Verilog-A compact model, compiled with OpenVAF to OSDI, and
run in ngspice, closing the physics → TCAD → SPICE loop with zero
commercial tools. Finally, we go below the semi-classical floor with a
self-consistent Schrödinger–Poisson solver (verified against the zeros of the
Airy function to 0.004 %, and with the Si/SiO₂ barrier treated as the finite
3.1 eV step it is, plus exchange–correlation), which places the
inversion-charge centroid ~2.1 nm from the interface — costing 4.1 % of the
inversion capacitance — and, fitted to a
measured 40 µm × 40 µm sky130 MOS varactor whose known area breaks the usual
thickness/area degeneracy, extracts a physical gate oxide of 3.60 nm against
the 4.15 nm the PDK's compact models carry, the difference being an
electrical thickness that silently contains the quantum centroid.
When Google open-sourced the SkyWater Sky130 PDK, an open question was posed by Tim Ansell in August 2021 on the DEVSIM issue tracker (devsim/devsim#74): could a fully open-source TCAD flow be brought up on Sky130 and shown to match the PDK's published data? The issue closed in April 2022 with zero comments. The ingredients have been public all along:
- the PDK's device menu with e-test targets and corner SPICE models,
- raw measured I-V data from a SkyWater test tile (google/skywater-pdk-sky130-raw-data, measured by CoolCAD Electronics; archived upstream April 2026 — clone a copy while you can), and, since 2026,
- a 77 K characterization dataset from UT-Arlington/Fermilab (repo, companion to arXiv:2604.21625).
What Sky130 does not publish is process information: doping profiles, halo/LDD implants, or any recipe data. The PDK maintainers have stated that reconstructing the process from public data is practically hopeless. This work therefore takes the device-physics route: geometry from the published stack, doping treated as fitted parameters, and measured silicon as the judge — which turns the fit itself into a measurement of the undisclosed quantities.
| Layer | Tool | Version |
|---|---|---|
| Device simulation | DEVSIM (pip install devsim + MKL) |
2.10.0 |
| From-scratch solvers | numpy / scipy | 2.5 / 1.18 |
| Compact-model compiler | OpenVAF-Reloaded | v24.0.1mob |
| Circuit simulation | ngspice (OSDI loader) | 46 |
| Linker for OpenVAF | llvm-mingw | 20260616 |
The full derivation, written to be readable, is in
01_pn_1d/EQUATIONS.md. In brief, all device
behavior here follows from three coupled PDEs — Poisson's equation and the
electron/hole continuity equations with drift-diffusion transport —
discretized with Scharfetter–Gummel fluxes [1] and Shockley–Read–Hall
recombination. On top of this baseline:
- Surface mobility (Sec. 3.2): µ = µ₀ / (1 + (|E⊥|/E₀)^α), with E⊥ the field component normal to the Si/SiO₂ interface, evaluated per mesh element and differentiated symbolically into the Newton Jacobian.
- Velocity saturation (Sec. 3.3): Caughey–Thomas form [2] with β = 2, using the channel-direction field component.
- Incomplete ionization (Sec. 3.4): at 77 K the acceptor charge in Poisson's equation becomes N_A⁻ = N_A / (1 + g_A·p/p₁A) with p₁A = N_V(T)·exp(−ΔE_A/kT), ΔE_A = 45 meV, g_A = 4 — i.e. dopant freeze-out (~3 % bulk ionization at equilibrium). Degenerate 10²⁰ cm⁻³ source/drain regions are kept fully ionized (Mott transition).
Three scalars are fitted per temperature — a gate workfunction shift (threshold position), a low-field mobility µ₀ (transconductance level), and the mobility-degradation exponent α (curve shape) — plus, where physics demands, one additional single-purpose parameter (channel doping from the body effect; halo doping from short-channel VT; saturation velocity from a single on-current point). Everything else reported is prediction. Fitted and predicted quantities are labeled in every results table below.
| Stage | Test | Result |
|---|---|---|
| 1 | from-scratch 1D pn diode vs closed-form theory | depletion width 2.9 % off analytic; ideality 1.015 (diffusion), 1.32 (SRH low bias — emergent, not programmed); current conservation exact |
| 2 | same diode, independent implementation (DEVSIM) | I-V ratio 1.000 over 5 decades; worst 1.9 %; equilibrium potential within 0.42 mV |
| 3 | 1D MOS capacitor (t_ox = 4.15 nm) vs textbook | V_T from ψ_s = 2φ_F matches the analytic formula to 0 mV; C-V regimes and HF C_min within 10 % |
Script: 04_sky130_nfet/mosfet_sky130.py,
mosfet_mobility.py.
Data: measured .mdm transfer curves (V_D = 0.1/1.8 V; V_B = 0/−0.9/−1.8 V).
| Quantity | Value | Status |
|---|---|---|
| V_T (max-g_m) | 0.516 V, matched to 0 mV | fitted (workfunction) |
| peak g_m | matched to 0.1 % | fitted (µ₀) |
| channel doping N_A | 3.94×10¹⁷ cm⁻³ | extracted from measured body effect (γ ∝ √N_A; initial 10¹⁷ guess predicted 153 mV shift vs 304 mV measured) |
| body effect ΔV_T(V_SB=1.8) | 310 mV vs 304 mV measured | predicted after N_A extraction |
| µ₀ / α / E₀ | 749 cm²/Vs / 1.96 / 5×10⁵ V/cm | fitted (shape) — note µ₀·E₀^α degeneracy: E₀ must be pinned and α fitted |
| linear I_D-V_G, on-region | mean ratio 0.982 (worst 0.969) | fitted region |
| saturation I_D-V_G (V_D=1.8) | mean ratio 0.938 (worst 0.920) | prediction — no saturation data used in the fit |
The constant-mobility baseline (no µ(E⊥) model) errs by +12 %/+22 % (mean/worst) in saturation; the field-dependent mobility halves-to-quarters every error:
Script: 04_sky130_nfet/mosfet_short.py.
All Sec. 3.1 parameters transferred unchanged; velocity saturation added.
The pure transfer is the experiment. With long-channel parameters only, the simulated V_T collapses to 0.297 V against 0.707 V measured — a −410 mV short-channel roll-off. Real Sky130 suppresses this with halo implants, which are not public; the model therefore detects their existence and measures their integrated strength: matching V_T requires an effective channel doping of 1.07×10¹⁸ cm⁻³ (secant fit, V_T matched to 2 mV) — 2.7× the long-channel value.
| Quantity | Value | Status |
|---|---|---|
| halo-proxy N_A | 1.07×10¹⁸ cm⁻³ | fitted (V_T only) |
| v_sat | 6.8×10⁶ cm/s | fitted (one point: I_on(1.8, 1.8) to 0.8 %) |
| DIBL (const-current V_T) | 121 mV vs 75 mV measured | prediction |
| saturation transfer (V_D=1.8) | mean ratio 1.045 | prediction |
| output curves V_G = 1.44 / 1.8 V | mean ratios 1.042 / 1.055 | prediction |
| output curve V_G = 1.08 V | 1.140 | prediction |
| linear region (V_D=0.1) | 1.615 → 1.045 after stage 4d (below) | was the documented limitation: ideal, perfectly conducting source/drain |
Stage 4d — series resistance (mosfet_rsd.py).
Two physical additions close the linear-region gap: (i) Masetti
doping-dependent mobility with literature parameters — normalized so the
calibrated channel is untouched, the 10²⁰ cm⁻³ n+ regions drop to 29 % of
channel mobility, no fitting — and (ii) a lumped contact resistance in
series with source and drain, solved self-consistently around the TCAD
device (internal V_S′ = I·R_s, V_D′ = V_D − I·R_d, damped fixed point).
One fitted value: R_c = 384 Ω·µm total, a textbook-plausible silicide
contact resistance. Results: linear region 1.615 → 1.045; saturation
1.015; all four output curves within 0.99–1.09. Once R_c stopped being
absorbed into other parameters, the fitted saturation velocity moved from
6.8×10⁶ to 8.3×10⁶ cm/s — toward the ideal silicon value of 10⁷: adding
the right physics made every other parameter more physical simultaneously.
Script: 04_sky130_nfet/mosfet_cryo.py.
Device: NMOS 100 µm × 100 µm, measured in liquid nitrogen at Fermilab
(dataset of arXiv:2604.21625).
Temperature enters from first principles: V_t = kT/q = 6.64 mV, n_i ≈ 2×10⁻²⁰ cm⁻³ (silicon is an insulator until gated), and dopant freeze-out — only ~3 % of channel acceptors are ionized at equilibrium (kT = 6.6 meV ≪ ΔE_A = 45 meV) — written into Poisson's charge density as a hole-dependent ionization term (Sec. 2.2). The same three-knob fit protocol as at 300 K is applied; everything else is prediction.
| Quantity | Value | Status |
|---|---|---|
| V_T | 0.607 V vs 0.606 V measured | fitted (workfunction; moved only +46 mV from its 300 K value — temperature physics carries the rest of the V_T shift) |
| µ₀(77 K) / µ₀(300 K) | 4990 / 749 ≈ 6.7× | fitted level; the ratio is a physical result (phonon freeze-out) |
| α at 77 K | rails at the 3.5 bound | the universal-mobility form structurally cannot reproduce the measured high-V_G droop — consistent with unmodeled Coulomb/trap scattering |
| saturation transfer (V_D=1.85) | mean ratio 1.27 | prediction |
| output V_G = 1.11 / 1.85 V | 1.33 / 2.0 | prediction; the V_G=1.85 excess is consistent with trap scattering + self-heating (~1.7 mW in LN₂), both absent from the isothermal model |
For scale: the dataset authors' own isothermal BSIM4 compact models, fitted directly to this data, report ~20 % relative RMS error — the physics-based predictions here land in the same class while additionally explaining the V_T shift and mobility change rather than refitting ~100 parameters.
The central finding — subthreshold-swing saturation, quantified. Boltzmann statistics require SS = ln(10)·n·kT/q. At 77 K the simulation — faithfully executing ideal physics, including the freeze-out-modified depletion capacitance — produces 23 mV/dec. The measured device: 71 mV/dec, barely steeper than its room-temperature value. This 3× gap is the well-known cryogenic SS-saturation anomaly, attributed to band-tail states and interface traps [3,4], and it is precisely the regime where standard drift-diffusion (and standard compact models) break. An ideal-physics simulation makes the anomaly measurable as a model residual:
(Numerical footnote: resolving a ~20 mV/dec slope requires a 5 mV gate sweep — the entire subthreshold region spans a few tens of mV and is invisible on a standard 50 mV grid.)
Ruling out series resistance (stage 4f,
mosfet_cryo_rsd.py). Could the
77 K on-state residuals just be the R_c = 384 Ω·µm extracted at 300 K? At
W = 100 µm that is 1.9 Ω per side, so the correction is perturbatively
small and first-order perturbation theory (ΔI/I = g_m·R_s + g_d·(R_s+R_d))
answers rigorously: the largest possible correction anywhere is 1.3 %,
against observed residuals of 27–100 %. Series resistance is ruled out;
the residuals belong to trap/Coulomb scattering and self-heating.
What causes the SS saturation? Model discrimination (stage 4g,
mosfet_cryo_tail.py). The
literature offers two mechanisms; both were tested.
- Hypothesis A — saturating acceptor interface-trap sheet (the classic D_it picture), implemented in Poisson's equation with occupation N_filled = N_t·(n/N_c)^(kT/W_t): falsified twice over. The simulated swing is set by the sheet density, not the tail energy — a 2.3× change in W_t leaves it at 185–197 mV/dec, 2.7× above measurement — and trapped charge is temperature-blind: the same sheet drags the already-validated 300 K curves by up to 43 % (q·N_t/C_ox ≈ 0.58 V of stretch-out at any temperature).
- Hypothesis B — band broadening (carrier statistics follow the tail, kT → W₀ [3]): survives every test this data allows. The extracted tail energy is W₀ = 20.5 meV (T₀ = 238 K) — squarely in the Urbach range; W₀ < kT(300 K) = 25.9 meV explains naturally why the 300 K device behaves ideally (the tail hides inside thermal broadening — exactly what the trap sheet could not do); and the ideal-physics curve with its voltage axis dilated by SS_meas/SS_ideal tracks the measured subthreshold to 42 mV RMS (0.59 decades) over 3.5 decades. The residual is genuine curvature — the real tail is not a perfect single-slope exponential.
Script: 05_compact/compact_fit.py,
model: 05_compact/ekv_lite.va.
A 4-parameter EKV-style model (single smooth expression across all operating regions) is fitted to the calibrated TCAD (4.6 % log-space RMS: VTO = 0.421 V, n = 1.218, KP = 261 µA/V²), compiled with OpenVAF to an OSDI shared library, and executed by ngspice 46 as a circuit device. The ngspice/OSDI implementation agrees with the reference implementation of the same equations to 0.7 % worst-case, and tracks the measured silicon to 8 % (linear) / 17 % (saturation) — inherited from the TCAD stage.
Scripts: 06_schrodinger/sp_solver.py,
validate_schrodinger.py,
sp_moscap.py,
sky130_cv.py; math in
06_schrodinger/EQUATIONS.md.
Everything above counts carriers with statistics applied to a local band edge. At a MOS surface that is wrong: the inversion layer is a ~1 nm well, motion perpendicular to the interface is quantized, and the electron gas is two-dimensional. Stage 6 replaces that assumption with an effective-mass Schrödinger equation (two electron ladders and two hole ladders of the (100) surface, exact Fermi–Dirac subband filling) coupled self-consistently to the same Poisson solver, on the same mesh, by a Trellakis predictor–corrector (~10 outer iterations, no damping).
The eigen solver is verified before use: an infinite square well to 0.003 %, a triangular well against the zeros of the Airy function to 0.004 %, its ground-state centroid to 0.002 %, and — the check that makes the comparison meaningful — a flat-band limit in which the sum over box subbands reproduces the classical bulk density to 0.07 %, still 0.11 % at a degenerate 8.8×10¹⁹ cm⁻³.
A methodological trap, recorded because it inflates the answer. The quantum branch must use Fermi–Dirac statistics — a 1.8 V inversion layer is degenerate. If the semi-classical branch is left on Boltzmann (as stages 1–5 legitimately are, where every density is non-degenerate), then "quantum vs semi-classical" silently measures the statistics as well as the confinement: Boltzmann over-counts by 6.7 % at 5×10¹⁸ cm⁻³ and by 88 % at 5×10¹⁹ cm⁻³, which is where a classical inversion layer peaks. Both branches here use identical statistics; making that correction cut the apparent quantum effect by about 20 %.
On the sky130 stack (t_ox = 4.15 nm, N_A = 10¹⁷ cm⁻³) the consequences at the 1.8 V rail are (quoted here for the bare hard-wall model, and revised by the refinement ladder further down): the charge centroid moves from 1.06 nm to 2.37 nm, the density vanishes at the interface where the classical one peaks (its maximum sits 1.35 nm deep), the quasi-static capacitance falls 6.9 %, and the threshold rises +53 mV at fixed inversion charge. The centroid is equivalent to an oxide 3.2 Å thicker — obtained two independent ways, from the centroids geometrically and by asking the semi-classical solver which t_ox best mimics the quantum C–V. Notably the shift is invisible to a band-bending threshold criterion (φ_s = 2φ_F moves −2 mV): it is a charge-criterion effect, because quantization changes how much charge a given bending buys, not the bending itself.
The same measured-silicon archive used for I–V also contains C–V. The target chosen here is a 40 µm × 40 µm MOS varactor (801 points, ±4 V) rather than the transistor split C–V beside it, and the reason decides what can be concluded: at 1600 µm² with a 160 µm perimeter, fringe and overlap are sub-percent, so the gate area is known — which breaks the t_ox/area degeneracy that makes C–V extraction elsewhere a shape comparison rather than a measurement. (The L = 0.15 µm transistor C–V was tried first and abandoned: a halo'd short channel is not a 1D object, and least-squares "solves" that by inflating overlap capacitance to 3 pF against the 0.79 pF the data directly measures.)
Dividing ε_ox by the measured 8.416 fF/µm² gives 4.10 nm — and that division is the semi-classical extraction, since it assumes the charge sits exactly at the interface. Fitting instead, with area pinned to geometry, poly depletion carried explicitly (it has its own Q² signature: the measured rail peaks at +2.4 V and rolls off 1.1 % out to +4 V) and both models scored on identical data points:
| extracted physical t_ox | rail residual RMS | fitted V_FB | |
|---|---|---|---|
| naive ε_ox/C_max (= the PDK's number) | 4.10 nm | — | — |
| semi-classical fit | 3.80 nm | 0.061 % | +0.62 V |
| Schrödinger–Poisson, hard wall | 3.60 nm | 0.124 % | −0.04 V |
| Schrödinger–Poisson, finite 3.1 eV barrier | 3.60 nm | 0.030 % | +0.23 V |
| Schrödinger–Poisson, + XC + oxide charge | 3.60 nm | 0.175 % | +0.59 V |
(measurement point-to-point noise: 0.016 %)
The extracted oxide is identical across every quantum variant — hard wall, finite barrier, exchange–correlation, oxide-resident charge — and across a 10× scan of well doping. The 3.60 nm does not depend on how refined the quantum model is, only on whether there is one at all. That invariance is the strongest part of the result; the rung-to-rung differences show up in the nuisance parameters and the residual, not in the answer.
The installed PDK's own numbers — BSIM toxe/toxm = 4.148 nm for
nfet_01v8 (with dtox = 0, i.e. it declares physical = electrical) and
cnwvc_tox = 4.165 nm for the varactor — agree with the naive extraction, as
they must: it is the same division. So the PDK's oxide thickness is an
electrical parameter that silently contains the inversion-layer centroid,
and this stage puts the physical oxide ~14 % thinner. The result is robust
where it matters: across a 10× scan of well doping the extracted t_ox does not
move at all (0.00 Å spread), against a 2.0 Å quantum correction.
The honest complication, and how it resolved. With a hard wall at the interface the quantum model does not fit better — its residual is 2× the semi-classical one, and both sit far above the 0.016 % measurement noise, so both are systematic model error. Confinement is not a free improvement in fit quality; what it changes is the physical parameter extracted.
That residual bowed, in the way an infinite barrier would if it overstated how fast the centroid shrinks with field, so the obvious suspect was the wall itself: a real Si/SiO₂ interface is a 3.1 eV step, not infinite, and the wavefunction leaks into the oxide. Rather than leave that as an explanation, it was run as an experiment — the Schrödinger domain extended across the oxide with a BenDaniel–Duke mass discontinuity (m_ox = 0.5 m₀), verified to reproduce the hard-wall limit as the offset is raised:
| hard wall | finite 3.1 eV barrier | |
|---|---|---|
| accumulation centroid ⟨z⟩ at matched bias | 2.483 nm | 2.316 nm |
| extracted physical t_ox | 3.60 nm | 3.60 nm |
| rail residual RMS | 0.124 % | 0.030 % |
The prediction survives: the gas moves in by 0.17 nm, the bow flattens, and the residual drops fourfold to become the best of all three models — twice the measurement noise rather than eight times it. And the headline is untouched: the extracted oxide does not move at all, so the 3.60 nm result does not rest on how the barrier is treated.
Two caveats stated rather than buried. This version leaves the ~0.1 % of the charge that sits inside the oxide out of Poisson, so it is not fully self-consistent. And the flatband voltage no longer ranks the models cleanly: V_FB should sit slightly below zero for n⁺ poly on an n-well, which the hard-wall fit gives (−0.04 V) and the semi-classical fit badly misses (+0.62 V), but the best-fitting barrier model lands at +0.23 V. The evidence there is mixed, and is reported as mixed.
06_schrodinger/refinements.py lifts the two remaining approximations and
measures every rung against the one below it, because corrections are neither
additive nor all worth having:
| rung | ⟨z⟩ (nm) | C/C_ox at 1.8 V | ΔC vs previous |
|---|---|---|---|
| semi-classical | 1.059 | 0.9205 | — |
| + quantization (hard wall) | 2.373 | 0.8567 | −6.93 % |
| + finite 3.1 eV barrier | 2.221 | 0.8655 | +1.03 % |
| + exchange–correlation (Hedin–Lundqvist LDA) | 2.071 | 0.8827 | +1.98 % |
| + oxide-resident charge in Poisson | 2.072 | 0.8828 | +0.02 % |
This revises a number published above. The 6.9 % capacitance loss quoted for the bare hard-wall model overstates the effect: the refinements claw ~3 % back, so the honest quantum correction on this stack is 4.1 %, with the centroid at 2.07 nm rather than 2.37 and the threshold shift +26 mV rather than +53. Exchange–correlation is the rung that mattered (−29 meV at 10¹⁹ cm⁻³, comparable to kT and to the subband spacing); the oxide-resident charge was worth 0.02 %, so the original approximation was justified — but that is now measured rather than assumed, which was the point.
A small physical bonus fell out of the LDA algebra: because Ry* ∝ m* and r_s ∝ m* cancel in Ry*/r_s, the exchange potential is set by the permittivity, not the effective mass — so the otherwise awkward choice of in-plane mass is worth only ~1 meV.
Remaining, untested: a proper 6-band k·p valence band — the two-ladder parabolic hole model is the crude part, so hole-side numbers stay indicative — and a 2D-specific exchange–correlation functional in place of the 3D LDA used here.
What this shows. All the pieces needed to answer devsim#74 exist and compose: an open device simulator, an open PDK with published measured silicon, and (new in 2026) an open cryogenic dataset. A three-parameter physics calibration on one measured curve yields few-percent predictions of independent curves at 300 K, transfers across a 167× change in channel length with one physically-motivated parameter each for halos and velocity saturation, and survives a 4× drop in thermal energy with one workfunction nudge and a mobility refit.
What was measured that Sky130 does not publish. Effective channel doping (3.9×10¹⁷ cm⁻³), effective halo strength (1.07×10¹⁸ cm⁻³), an effective surface-mobility law (µ₀ = 749 cm²/Vs, α ≈ 2 at E₀ = 5×10⁵ V/cm), a source/drain contact resistance (R_c ≈ 384 Ω·µm), a saturation velocity (8.3×10⁶ cm/s once R_c is separated out), and the 77 K mobility ratio (6.7×). These are effective parameters of a deliberately minimal model — not process values — but they are reproducible, and the protocol that produced them is in this repo.
What ideal physics cannot do (and why that is a result). The 77 K subthreshold measurement disagrees with Boltzmann statistics by 3×, and the 77 K high-field behavior overshoots any universal-mobility shape. Both residuals point at the same microscopic physics — band tails, interface traps, trapped-charge scattering — which is exactly the frontier the cryo-CMOS literature identifies [3,4]. A clean simulator turns those effects from "fitting difficulties" into quantified anomalies.
Limitations (all deliberate scope choices, all visible in the numbers): band-tail transport is identified (stage 4g) but not implemented inside the drift-diffusion statistics — the 77 K subthreshold is described, not self-consistently simulated; no self-heating (77 K high-power output curve); series resistance as a lumped element rather than silicide/contact geometry; no quantum confinement corrections (defensible at t_ox = 4.15 nm / 130 nm node, not below); uniform channel doping (halo modeled as an effective average); Boltzmann rather than Fermi–Dirac statistics.
Future work. Implementing band-tail (W₀) statistics inside the transport solver so the 77 K subthreshold is simulated rather than rescaled; a self-heating model for the cryogenic high-power regime; and — since the author's CORDIC-1 chip on this exact process arrives ~mid-2027 via TinyTapeout — closing the loop with self-measured silicon.
Runs natively on Windows 11 (no WSL); Linux/macOS should work with the same pip packages. Python ≥ 3.10.
pip install devsim mkl numpy scipy matplotlib
DEVSIM + pip-MKL note: the 2026 mkl wheel ships mkl_rt.3.dll, newer than
DEVSIM 2.10 probes for. Set DEVSIM_MATH_LIBS=mkl_rt.3.dll and add the
Python Library\bin directory to PATH — or just rely on
common/bootstrap.py, which every script here uses.
Measured data (gitignored, ~17 MB total — clone into refs/):
git clone --depth 1 https://github.com/google/skywater-pdk-sky130-raw-data refs/skywater-pdk-sky130-raw-data
git clone --depth 1 https://github.com/UTA-Advanced-Detector-Technologies/Skywater-130nm-77K-Cryogenic-Models refs/uta-77k
Stage-5 toolchain (gitignored, portable zips into tools/ — only needed for
the compact-model stage): ngspice-46 (Spice64), OpenVAF-Reloaded
v24.0.1mob (mingw build), llvm-mingw (provides the linker; copy
x86_64-w64-mingw32-gcc.exe to gcc.exe and set MSYSTEM=MINGW64 —
see the commit history for the full spelunking log).
Run order = stage order; each script is standalone, prints its pass/fail checks, and writes its figure next to itself:
python 00_smoke/cap1d.py # seconds
python 01_pn_1d/pn1d.py # ~1 min
python 02_devsim_diode/diode_devsim.py
python 03_moscap/moscap.py
python 04_sky130_nfet/mosfet_sky130.py # ~10 min (fit loop)
python 04_sky130_nfet/mosfet_mobility.py # ~15 min
python 04_sky130_nfet/mosfet_short.py # ~15 min
python 04_sky130_nfet/mosfet_rsd.py # ~20 min (needs short_results.npz)
python 04_sky130_nfet/mosfet_cryo.py # ~20 min
python 04_sky130_nfet/mosfet_cryo_rsd.py # seconds (perturbative analysis)
python 04_sky130_nfet/mosfet_cryo_tail.py # ~15 min (needs cryo_results.npz)
python 05_compact/compact_fit.py # ~1 min
python 06_schrodinger/validate_schrodinger.py # ~1 min (no measured data needed)
python 06_schrodinger/sp_moscap.py # ~4 min
python 06_schrodinger/sky130_cv.py # ~20 min first run, then cached
Stage 6 caches its solved sweeps in 06_schrodinger/.sweep_cache.npz
(gitignored): the tox x doping scan is ~44 self-consistent quantum sweeps on
the first run and instant on every run after. Delete the file to force a
re-solve — and note the cache key carries a physics version tag, so changing
the carrier statistics invalidates it rather than silently reusing stale
curves.
Hard-won numerical lessons are recorded where they were learned (commit
messages + code comments); highlights: DEVSIM leaves the corrupted Newton
iterate in the device after a convergence failure (snapshot/restore around
every solve, or bias bisection retries from garbage); rectangular-mesh
contacts need thin dummy regions behind the contact line; Newton RelError
plateaus at ~10⁻⁹ when carrier densities span 20 decades (tolerance ladder);
ramp the drain with the channel ON; ramp doping as a live parameter instead
of rebuilding; ngspice 46 silently ignores parenthesized OSDI .model
parameters; and in stage 6, a quantum-vs-classical comparison is only
meaningful if both branches share carrier statistics, band-edge densities and
mesh — otherwise it measures the bookkeeping (see §3.5).
[1] D. L. Scharfetter and H. K. Gummel, "Large-signal analysis of a silicon Read diode oscillator," IEEE Trans. Electron Devices 16(1), 64 (1969).
[2] D. M. Caughey and R. E. Thomas, "Carrier mobilities in silicon empirically related to doping and field," Proc. IEEE 55, 2192 (1967).
[3] H. Bohuslavskyi et al., "Cryogenic subthreshold swing saturation in FD-SOI MOSFETs described with band broadening," IEEE Electron Device Letters 40(5), 784 (2019).
[4] A. Beckers, F. Jazaeri, and C. Enz, "Cryogenic MOS transistor model," IEEE Trans. Electron Devices 65(9), 3617 (2018).
[5] F. Beall et al., "DC Cryogenic Modeling of Open-Source SkyWater 130nm MOSFETs at 77 K Using BSIM4," arXiv:2604.21625 (2026). Dataset: UTA-Advanced-Detector-Technologies/Skywater-130nm-77K-Cryogenic-Models.
[6] J. E. Sanchez, DEVSIM TCAD semiconductor device simulator, https://devsim.org (Apache-2.0).
[7] Google / SkyWater Technology, SKY130 open PDK and raw characterization data, https://github.com/google/skywater-pdk-sky130-raw-data (Apache-2.0, measurements by CoolCAD Electronics LLC).
[8] T. Ansell, "Building out fully open source TCAD for the SkyWater 130nm node," devsim/devsim#74 (2021).
All simulation code and result data in this repository; all measured silicon
data from the public repositories cited above ([7] is archived upstream —
mirror early). Code licensed Apache-2.0 (see LICENSE). To cite this work,
use the CITATION.cff file (GitHub → "Cite this repository").
Acknowledgments: measured data by CoolCAD Electronics/SkyWater/Google and by the UT-Arlington Advanced Detector Technologies group at Fermilab; DEVSIM by Juan E. Sanchez; ViennaTools, OpenVAF and ngspice communities.