# Coma, trefoil and spherical aberration in the model-based SSB fit (2026-10-21)

> **Historical:** describes the code as of 2026-10-21. The `higher_order=` subset fit it describes (a finite-difference Newton on the CPU with its coherence and precision guard) and its CUDA per-term agreement kernel were removed on 2026-10-25: the default fifth-order fit now fits coma, trefoil, C30 and every term through fifth order on exact gradients, with the half-band check, and `higher_order=` prints one line. The MPS per-term kernel remains as the fifth-order objective's reference. See the [current fit methods](../../../api/ssb.md).

## Question

The joint fit (`find_aberrations()`) maximizes the least-squares agreement of the
thick-sample SSB model with G(q, k) over C10, C12, phi12, the sample tilt and a
depth spread. The model's probe phase can carry any Krivanek term. Can the same
objective also recover coma (C21), trefoil (C23) and spherical aberration (C30),
jointly with the tilt and depth spread, accurately and fast? Coma and a sample
tilt both move bright-field images, so the first question is whether they can be
told apart.

The probe phase of the first sideband must be evaluated at k − q for the odd
orders; every backend had it at q − k, which `2026-10-21-ssb-odd-aberration-sign`
fixes and verifies. With the old sign, the model's agreement at the true
coefficients was 38% below the corrected model's, and the fit found no coma.

## Change

`find_aberrations(higher_order=("C21", "C23", "C30"))` (any subset; MPS and
CUDA) runs the usual fit: search, both 180-degree branches, optional rotation
refinement. Then one trust-region Newton polish runs from that fit. It moves
C10, C12 (as a vector), the tilt, the depth spread and the higher orders (each as
a Cartesian vector, C30 as a signed number), starting with the higher orders
at 0 (`direct_ptycho/higher_order_fit.py`). The objective is the joint fit's. A per-term
kernel evaluates it, computing the probe phase from the Krivanek polynomial
(`direct_ptycho/mps/higher_order_agreement.py`, `direct_ptycho/cuda/higher_order_agreement.py`); the
joint fit's angle-sum kernel needs a quadratic χ. `tilt=False` holds the tilt
and depth spread at 0.

The agreement evaluations and reconstruction run on the selected GPU backend.
The optional higher-order polish uses NumPy on the CPU for its Newton steps,
finite differences and Hessian calculations. It is not a GPU-resident optimizer.
The result records `host_numerical_search=True`; its `higher_order` record
identifies `optimizer_backend="numpy"` and the `agreement_backend`.
This optional path does not change the default fit optimizer.

The result's `aberrations` carry `C21_nm`, `phi21_deg`, `C23_nm`, `phi23_deg` and
`C30_nm`, so `reconstruct(result)` applies them. `result.optimizer["higher_order"]`
records the agreement before and after, the evaluations, the time, and the
guard below: `trusted` and its numbers in `guard`.

| kernel check | MPS | CUDA |
| --- | --- | --- |
| higher orders 0, against the joint fit's kernel | 1.1e-6 relative | 1.9e-6 |
| C21, C23, C30, against a float64 full-plane model (k − q) | 2.6e-7 (q − k: 6–40% off) | 2.6e-7 (q − k: 5–38% off) |
| ms per evaluation, 3,900 BF pixels, batch of 78 | 0.11 | 0.31 |

The CUDA timing was taken on a shared card. CUDA's 128 × 128 FFT kernels have no
higher-order path; a preview with higher orders and no depth spread now
reconstructs there through the depth-aware reference path at zero depth, which
is the same model and estimator.

## Evidence

Multislice simulations of an SrTiO3-type perovskite [001] supercell (6 × 6
cells, 23.4 Å; Gaussian atomic potentials, one frozen-phonon configuration),
300 kV, 30 mrad, 64 × 64 scan at 0.366 Å, 0.84 mrad detector pixels. Truth:
C10 −5 nm (mid-depth), C12 2 nm at 30°, C21 50 nm at 60°, C23 40 nm at −20°,
C30 3 µm, tilt (3, −4) mrad. The potential was scaled to separate the model's
identifiability from its validity: ×0.02 or ×0.1 (weak phase, where the SSB
model is exact), ×0.25 (a light-element crystal) and ×1 (SrTiO3).

**Identifiability.** At the fitted maximum, the correlations from the inverse
Hessian of the agreement (11 parameters) are:

| | weak phase, 15.6 nm, tilted | light, 8 nm, tilted | SrTiO3, 1 nm |
| --- | --- | --- | --- |
| coma with tilt | ≤ 0.07 | ≤ 0.23 | ≤ 0.24 (tilt not determined at 1 nm) |
| coma with C10, C12 | ≤ 0.26 | ≤ 0.15 | ≤ 0.05 |
| trefoil with anything | ≤ 0.19 | ≤ 0.18 | ≤ 0.06 |
| C30 with C10 | −0.85 | −0.86 | −0.93 |

The tilt enters the model only through the real depth weights, coma through the
phase of the transfer function, so they separate. Direct controls on the weak
15.6 nm crystal: data with the tilt and no aberrations above C12 fitted to a
coma of 1.3 nm, and holding the tilt at 0 left the agreement 16% below its
maximum with coma 2.7 nm (coma does not absorb the tilt); data with coma and no
tilt gave a tilt of (0.0, 0.0) mrad with or without coma in the fit (the tilt
does not absorb coma). C30 trades against C10 (balanced spherical aberration)
but is still determined.

**Accuracy, noise-free** (fitted magnitude; angles within 2° unless noted):

| crystal | C21 (50 nm) | C23 (40 nm) | C30 (3 µm) | tilt (3, −4) |
| --- | --- | --- | --- | --- |
| weak phase, one slice | 49.9 | 40.0 | 2.98 | not determined |
| weak phase, 15.6 nm | 50.3 | 39.9 | 3.08 | (3.0, −4.0) |
| light, 2 nm | 49.1 | 40.1 | 2.74 | (2.6, −2.3) |
| light, 8 nm | 46.9 | 39.2 | 2.81 | (3.1, −3.8) |
| light, 15.6 nm | 44.8 (15° off) | 38.5 | 1.04 | (3.0, −3.9) |
| SrTiO3, 0.4 nm | 48.7 | 39.9 | 2.80 | not determined |
| SrTiO3, 1 nm | 46.5 | 39.4 | 2.31 | not determined |
| SrTiO3, 2 nm | 39.8 | 37.5 | 1.80 | not determined |
| SrTiO3, 4 nm | 23.0 (9° off) | 32.4 | 0.45 | (−0.4, −1.7) |
| SrTiO3, 15.6 nm | 64.8 (25° off) | 19.2 | −3.45 | (3.4, 0.6) |

Where the crystal scatters weakly the fitted terms are the probe's to a few
percent. In strongly scattering crystals they are not: an SrTiO3 crystal of
15.6 nm with no aberrations above C12 fitted to coma 55 nm and C30 −7 µm, while
the agreement rose 11%. The terms then fit dynamical scattering the model does
not contain. A coma of 100 nm, trefoil of 100 nm and C30 of 10 µm (light, 4 nm)
were recovered as 98 nm, 100 nm and 9.7 µm from a joint fit whose tilt was
8 mrad off; 10 nm, 10 nm and 1 µm as 9.8 nm, 10.2 nm and 0.95 µm.

Fitting the higher orders also corrects the low orders: with the truth above,
the joint fit alone was 0.6 to 1.4 nm off in C10 (C30 balanced by defocus) and
1.2 mrad off in tilt (light, 8 nm); with them, 0.2 nm and 0.2 mrad.

**Noise.** Poisson noise at 1e5, 1e4 and 1e3 e/Ų (three seeds each) on the
64 × 64 scans scattered the fitted coma by 0.1 to 2 nm (standard deviation) and
C30 by 0.01 to 0.2 µm; a 10 nm coma's azimuth, 5° off without noise, was
18 ± 7° off at 1e3 e/Ų. For these crystals the bias above, not the dose,
limits the accuracy. On the weak-phase 15.6 nm crystal (little
signal) 1e4 e/Ų scattered the coma by 3.5 nm and 10°.

**Cost.** One polish: 150 to 1,400 evaluations in 4 to 21 batched passes. On the
64 × 64 scans (3,900 bright-field pixels) 0.04 to 0.4 s, on top of a joint fit
of 0.4 to 2 s. On another 512 × 512 acquisition (8,807 bright-field pixels) the
polish took 7.2 s (437 evaluations, 8 passes, 16 ms each) after a two-branch
joint fit of 10.3 s; with the rotation given the whole fit went from 7.3 s to
13.0 s (C21, C23, C30), 12.0 s (C21, C23) and 8.5 s (C30). On 256 × 256 regions
it added 1.7 to 2.4 s. M5 Max.

## CUDA and MPS give the same fit

The same 128 × 128 scans (the 64 × 64 simulations above, tiled; one period, so
neither backend pads) were fitted on both backends with `tilt=True`:

| crystal | backend | C21 | C23 | C30 | tilt (mrad) | polish |
| --- | --- | --- | --- | --- | --- | --- |
| weak phase, one slice | MPS | 50.08 nm @ 60.00° | 40.09 nm | 2911.8 nm | (0.33, 0.07) | 0.39 s |
| | CUDA | 50.08 nm @ 60.00° | 40.09 nm | 2911.5 nm | (0.33, 0.06) | 0.28 s |
| weak phase, 15.6 nm | MPS | 50.30 nm @ 61.68° | 39.93 nm | 3080.1 nm | (2.98, −4.03) | 0.20 s |
| | CUDA | 50.30 nm @ 61.68° | 39.93 nm | 3080.0 nm | (2.98, −4.03) | 0.20 s |
| light, 8 nm | MPS | 46.87 nm @ 58.14° | 39.23 nm | 2822.3 nm | (3.12, −3.80) | 0.22 s |
| | CUDA | 46.87 nm @ 58.13° | 39.23 nm | 2825.1 nm | (3.12, −3.80) | 0.08 s |

The agreements at the two maxima agree to 1e-7 relative. With `tilt=False`
(weak lattice) C21 was 49.97 nm @ 60.00° on both backends. On the measured 512 × 512
acquisition (8,807 bright-field pixels, both 180-degree branches) the two
backends returned the same terms to 0.01 nm and 0.01° (C21 104.45 against
104.46 nm, C30 5276.6 nm) and the same tilt and depth spread to 1e-5. The CUDA
polish took 4.9 s (437 evaluations in 8 passes, as on MPS) and the whole fit 12.1 s, against
7.2 s without the higher orders; on MPS, 7.1–7.4 s, 17.7–20.9 s and 11.3 s.
CUDA ran on an RTX PRO 6000 Blackwell card that other jobs kept at 99%
utilisation; MPS ran on an M5 Max.

## When to trust the terms

A fitted term is the probe's only where the SSB model describes the data. Five
guards were tested on 280 fits of new 64 × 64 multislice simulations, each
noise-free and at 1e5, 1e4 and 1e3 e/Ų with three Poisson draws:

- **Thick nulls:** 12 crystals with no aberration above C12 in the truth.
  SrTiO3-type at 2, 4 (tilted and not), 8 (two settings), 15.6 (tilted and
  not) and 30 nm; the potential ×1.5 at 4 and 8 nm; ×0.25 at 15.6 and 30 nm.
- **Weak-phase cases:** 9 crystals, potential ×0.1 or ×0.02, 0.2–30 nm. They
  carry coma, trefoil and C30 of 50/40 nm and 3 µm, of 10/10 nm and 1 µm, of
  25/20 nm and −2 µm, coma alone (30 nm), and none.
- **Thin series:** the light-element and SrTiO3 crystals of the earlier study
  (2–8 nm light, 0.4–1 nm SrTiO3).

A fit counts as accurate when coma and trefoil are within max(5 nm, 15%)
(vector error) and C30 within max(0.5 µm, 25%) of the truth. Without a guard,
100 of the 120 null fits are wrong: the terms absorb dynamical scattering (coma
up to 200 nm, C30 to −15 µm). So are 30 of the 90 weak-phase fits, all at
1e3–1e4 e/Ų on the weakest specimens, where noise sets the error.

| guard | thick nulls passed | weak phase passed (wrong among them) | accurate weak passed |
| --- | --- | --- | --- |
| none | 120 / 120 | 90 / 90 (30) | 60 / 60 |
| agreement gain over the fit without higher orders ≥ 1% | 51 / 120 | 59 / 90 (24) | 35 / 60 |
| depth spread inside the depth model's validity range | 118 / 120 | 90 / 90 (30) | 60 / 60 |
| scatter over the three noise draws ≤ 2 nm (coma) | 84 / 120 | 36 / 90 (0) | 36 / 60 |
| coma and trefoil agree within 5 nm between two q bands | 13 / 120 | 22 / 90 (0) | 22 / 60 |
| coherence ≥ 0.9 | 0 / 120 | 58 / 90 (7) | 51 / 60 |
| **coherence ≥ 0.9 and edge-phase error ≤ 0.1 rad** | **0 / 120** | **45 / 90 (0)** | **45 / 60** |

The depth model's validity range is 50 nm for Si, 30 nm for SrTiO3 and 20 nm for
BaTiO3, taken as ×0.25, ×1 and ×1.5 here. Neither the agreement gain nor that
range separates the nulls: a thick null gains as much from the extra terms as a
small true coma does, and the higher orders fail far inside the depth model's
range (SrTiO3 from 4 nm). Scatter over noise draws and the two-band split do not
see a bias that every draw shares.

**Coherence** is the share of G's signal power over the fit band that the fitted
model explains:

ρ = (A − F_A) / (P − F_P), with P = Σ_q w_q Σ_k |G(q, k)|².

A is the agreement and A ≤ P. F_A and F_P are their noise floors; for counted
electrons the noise of G(q, k) has variance G(0, k). A weak-phase object gives
ρ ≈ 1, because its G is the model's. Dynamical scattering puts power where the
model cannot follow, and the extra terms then fit part of it. Noise-free,
the weak-phase crystals reached 0.96–1.00 and the thin series 0.90–0.99. No
null whose fit was wrong exceeded 0.80: light 15.6 nm reached 0.80, SrTiO3 4 nm
0.77, 8 nm 0.58 and 15.6 nm 0.25. The 2 nm SrTiO3 null, whose fit returned no
spurious terms, reached 0.87.

**Precision.** The covariance of the fitted coordinates is 2σ²(−∇²A)⁻¹, taken
from the polish's own Hessian. A / (2σ²) is the profile log-likelihood for
Gaussian noise of variance σ² on G over the independent half of the q plane.
The guard asks that each term's standard error change the probe phase at the
aperture edge by ≤ 0.1 rad: (2π/λ) σ_C α^(n+1)/(n+1) ≤ 0.1, which at 30 mrad
is 3.5 nm for coma and trefoil and 0.15 µm for C30. Across 79 (crystal, dose)
groups, the coma's scatter over three draws divided by this standard error had
a median of 1.00 (quartiles 0.74–1.38). The very noisy weak cases are the
exception: their fits leave the quadratic region, the scatter exceeds the
prediction, and the guard rejects them anyway.

With both conditions, no wrong fit of the 280 passed. Weak-phase fits passed
9/9 noise-free, 21/24 accurate ones at 1e5 e/Ų, 15/19 at 1e4 and 0/8 at 1e3.
The thin series passed 51/70, all accurate. The guard is conservative: the 20
accurate null fits (SrTiO3 2 nm, and 4 nm untilted) fail it at ρ = 0.74–0.87.

On the measured 512 × 512 acquisition (a thick crystalline film) the polish is
statistically precise (coma ±0.07 nm), but the model explains 0.27 of the
signal, so the guard rejects it. That fits the earlier finding that the coma
there varies by ±50 nm between quadrants that share one probe.

`find_aberrations(higher_order=...)` records `coherence`, `standard_error_nm`,
`edge_phase_error_rad` and `trusted` in `result.optimizer["higher_order"]["guard"]`.
An untrusted fit prints one line saying why and keeps the terms.
Input that is not integer counts gets no noise floor (F = 0), which lowers the
coherence of noisy data, and no precision estimate.

## Conclusion

Coma and trefoil are identifiable and separable from the sample tilt; C30 is
identifiable but strongly correlated with C10. The fit is accurate where the
SSB model holds (thin or weakly scattering specimens) and returns non-probe
terms in thick, strongly scattering crystals. It is therefore an option for
calibrating the probe on suitable specimens, not a default. The coherence and
precision guard tells the two apart without the truth: it passed no wrong fit
in these simulations, and it passes most weak-phase fits at 1e4–1e5 e/Ų.