The solver

This work develops numerical methods for the quasiclassical theory of superconductivity and superfluid ³He, 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. This is a cited theorem of the classical matrix Riccati theory (Reid, 1972), not a numerical observation; its application to these coefficients, with the derivations extending it, is the repository's standing record at plans/Matrix_Riccati_Ball_Bound.md. 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 (the closest published peer runs on CPU-cluster compute and claims no GPU execution; the closest GPU peer is singlet/spin-degenerate on its own scope statement).
• 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.

Verification

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 battery establishes — and, with equal care, what it does not cover — is set out in VERIFICATION.md.

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, in the singlet limit, 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⁻⁴%.
• Independent RK-45 Integrator: An independent adaptive Runge–Kutta integration of the nonlinear Riccati trajectory reproduces the finite-element forward field to a relative L₂ error of 2.91 × 10⁻⁷, the precision floor of the archived field export — an upper bound on the disagreement rather than a resolved solver error. Under a registered protocol with the order parameter prescribed and held fixed across meshes, the discontinuous-Galerkin transport solve measures convergence order p = 2.00 against an analytic-profile RK–45 reference.
• Linearized Control: A control solve with the quadratic term dropped misses the bulk amplitude by more than an order of magnitude — a verification instrument built to fail, distinct from the published one-factor-lagged scheme, which converges.
• Free Energy Evaluation: On a uniform state, GPU evaluation of the Luttinger–Ward functional matches the closed-form analytic result to three parts in 10⁸.

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, with the order parameter carried at full linear order within each cell and both Kronecker product orderings in the Newton step. It runs the full self-consistent loop to convergence on device at production resolution, its boundary excess at the continuous-Galerkin floor, with the continuous-Galerkin prescribed-profile study as the baseline beside its measured convergence order.

One finding from that auditing is open, and it is a property of the discretization rather than of any boundary condition. Along the trajectories directed normal to the walls, the computed coherence amplitude exceeds its physical bound — the closed unit disk — at the nodes where the pairing interaction steps. The excess is present when no reflecting term is assembled at all, so it is not a consequence of one: it predates the reflecting-wall work and is present in every inhomogeneous run this solver has performed, unrecorded until now only because the exported field slice had always been taken along a single trajectory direction, and the work on reflecting walls is what produced the diagnostic that exposed it. Its size tracks the element aspect ratio; a prediction drawn from that dependence was pre-registered, tested on a further mesh, and failed, while the dependence itself survived. The mechanism is not established.

Direction

The central challenge in thin films and slabs (D ~ 10 ξ₀; Vorontsov, Phil. Trans. R. Soc. A 376, 20150144) is that the coherence length ξ₀ represents both the scale of variation and the resolution limit of quasiclassical theory itself. Experimental nanofluidic regimes (D/ξ₀ = 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. The program's definition of complete is exactly that target — reproducing the confined chiral-³He results of Heikkinen et al., phase ranking by free energy included. The plan of record sequences the work toward it in six phases (plans/Solver_Completion_Plan.md §3), two of which sit off the ³He path itself: a benchmark against published d-wave results, and a temperature-sweep driver for the thermodynamics. The four capabilities the target needs, in a strict order that avoids confounding physical and numerical effects:

1. Carrying specular boundary conditions through the discontinuous Galerkin formulation — the perfectly specular limit. This step is built: the reflected-flux term, the perfectly reflecting limit of the surface-scattering condition, imposed on the inflow faces where the trace is a degree of freedom of the scheme, with a lagged once-per-iteration update supplying the reflected trace and held outside the self-consistency mixer's own history. It is gated at 218 checks over nine mesh-and-angle configurations, zero failures, host-side and again against a build for the hardware, with an identity map and a wrong map kept in the gate as permanent controls so that what the checks discriminate is the pairing itself. With the walls switched on, runs converge in the same iteration counts as their walls-off baselines. The target experiment's walls are near-specular rather than perfect (specularity S > 0.97, inferred from the small T_c suppression), and reproducing that paper's own gap-suppression calculation, computed at S = 0.98, would need a partial-specularity boundary condition, S < 1, which this step does not supply and is recorded as a named capability gap.
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.
3. Exercising physics in the non-commuting triplet sector.
4. Evaluating competing free-energy candidates to rank order-parameter configurations.

One design decision on that path is settled: strong-coupling corrections — required to split the chiral A phase from its weak-coupling-degenerate planar partner — enter through the free energy, added on the converged weak-coupling state, with the transport equations left unmodified; the temperature dependence of the published correction coefficients, derived near T_c, is carried as a bracketed input assumption rather than a single adopted rule. None of it is implemented: this is the recorded design for steps 3–4, grounded in the strong-coupling literature (Rainer & Serene, Phys. Rev. B 13, 4745; the Supplemental Material of Heikkinen et al.).

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.

Correspondence welcome — elizabeth@bosonuum.com