Research

This work develops numerical methods for the quasiclassical theory of superconductivity and superfluid 3He, with a focus on confined geometries (films and slabs) where surface effects reshape the entire sample. The core formulation uses the Eilenberger–Riccati transport equations, implemented as a finite-element solver, self-consistent on the Matsubara axis, executed on GPU.

It started as a physics problem and turned into a numerics problem: advection-dominated transport, characteristic sweeps, and stabilization at gradient layers—the same machinery used in aerospace, pointed at a superconducting coherence amplitude instead of a fluid velocity.

Method

Quasiclassical theory reduces the Gor'kov equations to transport equations along classical trajectories on the Fermi surface by integrating out variation on the scale of the Fermi wavelength while retaining variation on the scale of the coherence length.

The Riccati Parametrization
The Riccati parametrization rewrites these transport equations in terms of coherence amplitudes:

  • Normalization: The parametrization satisfies the normalization condition identically by construction, preventing discretization error from violating it.
  • Boundedness: Boundedness is a property preserved by the nonlinear Riccati flow—the trajectory forms a Möbius map carrying the unit ball into itself. A solution that drifts outside this bound is recognizably wrong without reference to external data.

Spin-Matrix Closure
Absorbing the particle–hole index leaves a 2 × 2 spin matrix carrying one singlet and three triplet components—all four required for multi-component order parameters. In this non-commuting sector, components couple and the Newton step acquires a Kronecker-structured Jacobian.

Contributions Over Published Record
Building on the scalar and spin-matrix discontinuous Galerkin formulation of Seja and Löfwander (Phys. Rev. B 106, 144511; Phys. Rev. B 110, 064502), this implementation adds:

  • GPU Acceleration: Execution of the full matrix closure on GPU (where published peers run on CPU clusters or remain singlet-restricted).
  • Full Newton Linearization: A full Newton step through a Kronecker-structured Jacobian carrying both product orderings with no factor lagged—and verified against central finite differences of the residual it linearizes.
  • Bounded-Flow Auditing: Production-run artifacts are audited against the amplitude's unit bound—an audit, not a clamp: the solver itself is untouched.
  • External Verification: Independent cross-validation against non-internal reference solutions.

Verification & Open-Source Practice

This code carries the GNU LGPL v3 — COPYING and COPYING.LESSER ship in the repository — ahead of its open-source release. Every term is separately gated, every verification number regenerates from a clean checkout, and the record of what broke ships alongside what works. What the verification battery establishes — and, stated with equal care, what it does not cover — is written for a reader who was not in the room in VERIFICATION.md at the repository root.

In practice, theory and execution are aligned by auditing the invariant directly on real runs via Reports/audit_gamma_unit_disk.py. Verification is conducted against references external to the solver code:

  • Exact Analytic Root: In a homogeneous s-wave reservoir, the solver converges the bulk gap onto the exact discrete BCS root to within 1.3 × 10−4%.
  • Independent RK-45 Integrator: An independent adaptive Runge–Kutta integration of the nonlinear Riccati trajectory reproduces the finite-element forward field to a relative L2 error of 2.91 × 10−7 — the precision floor of the archived six-digit field export, an upper bound on the disagreement rather than a resolved solver error. A four-rung mesh ladder on the discontinuous-Galerkin path (2026-08-21) tightened the bound at full export precision: agreement within 6.5 × 10−9 at every resolution — at, and mostly below, the reference integrator's own tolerance (10−8), so most of the measured disagreement is the reference's error, not the solver's. These are bounds, and were deliberately reported as bounds: that measurement was floor-limited by the reference itself, and the repository records why. A subsequent registered protocol (2026-08-23) converted the bound into a measurement: with the order parameter prescribed and held fixed across meshes, the discontinuous-Galerkin transport solve measures convergence order p = 2.00 (full precision 1.9991456569968427) against an analytic-profile RK–45 reference at 10−12 relative tolerance, the measurement's floors recorded beside every reading in the repository's run record.
  • Linearized Control: A control solve with the quadratic term dropped misses the bulk amplitude by more than an order of magnitude.
  • Free Energy Evaluation: On a uniform state, GPU evaluation of the Luttinger–Ward functional matches the closed-form analytic result to three parts in 108.

Results to date lie in the singlet sector: no triplet state is constructed, and the triplet physics remains unexercised—the non-commuting algebra is verified structurally, not physically. A discontinuous-Galerkin formulation sharing the same closure is implemented and separately verified term by term — since 2026-08-19 with the order parameter carried at full linear order within each cell and both Kronecker product orderings in the Newton step, the device-side sizing and layout gates passing their first GPU builds at 24 and 16 checks, zero failures — and has run the full self-consistent loop to convergence on device at production resolution (2026-08-21), its boundary excess at the continuous-Galerkin floor. The convergence-order comparison against the continuous baseline remains open, and is reported as open.

Direction

The central challenge in thin films and slabs (D ~ 10 ξ0; Vorontsov, Phil. Trans. R. Soc. A 376, 20150144) is that the coherence length ξ0 represents both the scale of variation and the resolution limit of quasiclassical theory itself. Experimental nanofluidic regimes (D/ξ0 = 1; Heikkinen et al., Phys. Rev. Lett. 134, 136001) leave no region untouched by boundary conditions.

Interpreting these experiments requires theory that is quasiclassical, matrix-valued, and free-energy-resolved. Development beyond the open-source release — re-scoped in August 2026 to follow the release rather than gate it, and kept as the map — follows a strict four-step sequence to avoid confounding physical and numerical effects:

  1. Carrying specular boundary conditions through the discontinuous Galerkin formulation.
  2. Expanding from single-cell-deep strips to fully two-dimensional domains. Readiness is established: a two-row mesh runs on device as a pure mesh input with no code change, its two cell rows identical to machine precision (2026-08-21).
  3. Exercising physics in the non-commuting triplet sector.
  4. Evaluating competing free-energy candidates to rank order-parameter configurations.

Background

B.S. Aerospace Engineering and Mechanics, University of Minnesota; thirty years of independent consulting in computational modeling and systems analysis across thermal, fluid, and acoustic systems. Graduate Research Associate, Montana State University (2011–2012), implementing CUDA parallel architectures for fluid-flow simulation.

Open-source release under LGPL v3 — the license already ships in the repository. Correspondence welcome — elizabeth@elizabethburnim.com