Research

I build numerical solvers for the quasiclassical theory of superconductivity — the Eilenberger–Riccati transport equations that describe superconductors and superfluid ³He at scales between the atomic and the macroscopic.

The problem

A superconductor is described by an order parameter, Δ, that measures the strength and phase of electron pairing at every point in the material. In an infinite, uniform sample it is a constant and the physics is textbook. Near a surface it is not.

Within a few coherence lengths of a boundary, quasiparticles scatter in a way that converts electrons into holes and back — Andreev reflection. The states this creates sit inside the superconducting energy gap, and they restructure the low-energy spectrum of the material. In a sample only a few coherence lengths across, that restructuring is not a surface correction. It is the whole sample. The free-energy landscape changes shape, and phases become stable that cannot exist in bulk.

This is where the interesting physics is, and it is where the equations are hardest to solve.

What I am building

The quasiclassical theory reduces the full microscopic problem to transport equations along classical trajectories across the Fermi surface. The Riccati parametrization rewrites those equations in terms of coherence amplitudes that satisfy a normalization constraint identically, so the constraint cannot be violated by numerical error. The price is a nonlinear equation.

Seja and Löfwander (Phys. Rev. B 106, 144511) solve the scalar case with a discontinuous Galerkin finite-element method, self-consistently and for general device geometries. I am extending that formulation in two directions:

  1. Matrix-valued coherence amplitudes. The scalar case describes conventional spin-singlet pairing. The Riccati parametrization absorbs the particle–hole index, so what survives on the coherence amplitude is a 2×2 matrix in spin — one singlet and three triplet components. Carrying all four is a prerequisite for the multi-component order parameters of superfluid ³He and unconventional superconductors, and it is a generalization Seja and Löfwander explicitly defer.

  2. GPU execution. The problem is embarrassingly parallel across Fermi-surface direction and Matsubara frequency, and the arithmetic is double precision throughout.

The results below come from a continuous-Galerkin implementation of the same closure, stabilized by SUPG. The discontinuous formulation is written and in validation; a convergence comparison between the two is deferred to separate work.

Verification

A solver that converges is not a solver that is correct. This one is checked against references outside itself:

  • The converged bulk gap matches the exact analytic BCS root to 1.3 × 10⁻⁶.
  • The coherence field agrees with an independent Runge–Kutta integration, sharing no code with the solver, to a relative L₂ error of 2.9 × 10⁻⁷.
  • A deliberately failing control — the same comparison with the nonlinear term removed — establishes that the test is sensitive to the physics it claims to verify, rather than passing because both codes make the same simplification.

Applied to a superconductor–normal-metal–superconductor junction, the solver reproduces the canonical proximity profile: the order parameter vanishes in the normal region and heals back to the bulk gap over a structured shoulder, while the induced pair amplitude leaks continuously across the interface.

The matrix machinery is implemented and assumes no commutativity between the order parameter and the coherence amplitude. The results reported here lie in the singlet sector, where the two commute and the physics reduces to the scalar case — which is precisely why an exact analytic reference is available to check against, and why these results do not yet constitute a test of the triplet path.

One result worth recording is a failure. An early version advanced the forward and backward coherence fields under different closures, breaking a symmetry of the continuum theory at the level of the discretization. The gap settled into a smooth, plausible, wrong answer — 13.3% high, flat across the whole domain to five decimal places, with no boundary-layer structure to betray it.

What caught it was not a gate. It was reading the gap value against the exact analytic BCS root during runs launched to check something unrelated. No automated check for that class of error existed at the time; one was written four days later, after the fix.

The lesson is narrower than "the failure was silent," and more useful. The residual was not silent — it sat several hundred times above tolerance and the run never converged. What was silent was the physics: a spatially uniform, stable-looking order parameter, the shape the answer was supposed to have, sitting on top of an iteration that had not finished. An output can look settled while the diagnostics say otherwise, and the plausible shape is exactly what makes the diagnostics easy to discount. Symmetries of the theory have to be respected by the solver, not merely by the equations it discretizes — and the thing that reveals it when they are not is a reference the code cannot influence.

Where it is going

The target is confinement. Anton Vorontsov's review of Andreev bound states in films and slabs (Phil. Trans. R. Soc. A 376, 20150144) sets out what happens when a sample is squeezed to D ∼ 10 ξ₀: transitions that are absent in bulk, phases that break time-reversal or translational symmetry, and — in superfluid ³He — a competition between the B, planar and stripe phases decided by the structure of a multi-component order parameter rather than by its magnitude.

Recent nanofluidic experiments have reached D/ξ₀ = 1 (Heikkinen et al., Phys. Rev. Lett. 134, 136001). The theory needed to interpret them is quasiclassical, matrix-valued, and free-energy-resolved. The first two are what I am building.

Background

I came to this from aerospace engineering. The numerical machinery transfers more directly than the fields suggest: advection-dominated transport, characteristic sweeps, stabilization at gradient layers, and the discipline of verifying a code against something that is not itself. A superconducting coherence amplitude streaming along a Fermi trajectory obeys the same class of equation as a scalar advected by a flow.

Earlier work: parallel computing architectures for fluid simulation (Montana State University), computational modeling and test-system design in industry, and wind-tunnel aerodynamics.

Manuscript in preparation. Correspondence welcome — elizabeth@elizabethburnim.com