TDGL Simulation

tdgl3d is a fully three-dimensional simulator for superconducting devices and circuits, built on the time-dependent Ginzburg–Landau (TDGL) equations. It resolves what circuit models leave out: where vortices enter, how they move under a bias current, where they dissipate, and how much the film heats. Stacked metal–insulator–metal layers, holes, notches and the vacuum around the conductor are all on one 3-D grid, so the third dimension is part of the physics rather than an extrusion of a 2-D film.

The order parameter ψ and the gauge field live on a structured finite-difference grid, with A stored as link variables (Peierls phases) on the edges so that gauge invariance and flux quantisation hold exactly on the lattice rather than approximately.

Run it live in your browser → A 3-D NbN detector block on SiO₂/Si with bias current, field, heat, photon absorption and a ramped I–V curve, solved on your GPU (more below).

Everything is dimensionless: lengths in coherence lengths ξ, fields in $B_{c2}$. Three numbers fix the normalisation — $\Phi_0 = 2\pi$, $\lambda = \kappa$ (in ξ), and $H_{c2} = 1$. Devices drawn in SI units are converted with ξ taken at the temperature of interest.

What it simulates

• Vortex dynamics in 3-D — entry, pinning on holes, lattice formation and flux flow, with vortex lines free to bend and tilt through the thickness.
• Heterostructures — S/I/S stacks, inclusions, etched (vacuum) holes and regions of locally suppressed $T_c$, coupled only through the field.
• Transport current — normal-metal terminals inject current, with the electric potential solved at every step. Either the film is thin and does not screen (pyTDGL’s limit), or the current’s own magnetic field is carried self-consistently.
• Self-heating — a heat equation driven by the TDGL dissipation feeds back through the local temperature, so a hotspot can switch a wire and latch it normal.
• Single-photon detection — an absorbed photon is a hot spot in that heat equation, and the order parameter relaxes at the ratio u of a real superconductor, so detection thresholds and latencies of nanowire detectors come out of the same model.
• Heater-switched devices — resistive heaters, insulators that pass heat but no current, per-material conductivity and thermal properties, and a readout circuit (bias, load, series inductance) at the device’s terminals: enough to run a heater nanocryotron (hTron) end to end.
• Circuit parameters — kinetic and geometric inductance per unit length and small-signal AC impedance from a TDGL run, a static London solver for the port inductance matrix of 3-D wiring, and a 2-D electrostatic solve for capacitance, giving $Z_0$ and phase velocity.
• A service around it — every capability is described by a versioned project schema and runs as a job behind a FastAPI server, so a UI or an AI assistant drives the same simulations as a script.

In three dimensions

Runs where the thickness and the vacuum around the conductor change the answer. Each starts as a rendered image; Explore in 3-D loads the scene itself, which can be rotated, panned and zoomed, and a legend entry clicked to hide a layer.

A moat grid in a field-cooled S/I/S ground plane

A 36 × 36 ξ ground plane — 3 ξ metal, 1 ξ oxide, 3 ξ metal — with nine 4 ξ moats etched through it, cooled through T_c in B_z = 0.25 H_c2, about 52 flux quanta over the plane. The two metal layers are coupled only through the field and nucleate from independent random starts, yet 16 of the bottom layer’s 18 film vortices end up stacked on a top-layer vortex within 1.6 ξ. Each layer carries 39 quanta in all: the bottom keeps two more in the film, and the top holds those two in its moats instead.

A 3×3 array of 4 µm holes

The coherence length fixes the grid spacing, but fabrication fixes the device, so a real hole array is a large simulation: nine 4 µm holes on an 8 µm pitch with an 8 µm buffer is 36 µm across, and at ξ = 150 nm that is 240 × 240 × 9 — 457 k nodes. At ξ = 50 nm it is 15 M nodes in 12 GB.

The obvious way to get flux into the holes does not work, and that is the interesting part. Ramping the field up, nothing enters until 3.15 mT — and just above that, hundreds of vortices enter at once. Held at 3.6 mT, 567 of them pack the buffer into a triangular lattice while the array stays fully Meissner-screened behind them: the flux front stalls at the array perimeter and never reaches a hole. There is no applied field at which this film holds one or two vortices in equilibrium; it holds none, or it holds hundreds.

Field-cooling does what the experiment does. ψ grows from near zero with the field already on, so flux is trapped where it is rather than having to cross 8 µm of screening metal; the field then drops below the entry threshold and the state settles. Cooled at 4.0 mT and held at 2.0 mT, after 400 τ_GL: 3 vortices in the metal between the holes, 7 flux quanta trapped across the nine, with the rest of the flux in the buffer lattice. An independent run of the same protocol at another noise seed gives 2 and 6, with the census flat from t = 293 to t = 400 — so this is a settled state, not one realisation of the noise.

The holes hold about one quantum each rather than the $B A/\Phi_0 \approx 19$ the applied field would suggest, and the $B_z$ map says why: 8 µm of buffer at λ = 300 nm leaves the array interior nearly field-free, so there is no local field there to support more. Per-hole occupancy is set by how well the surround screens, not by the applied field — which makes the buffer width and κ the levers.

A vortex in the film and a fluxoid in a hole are counted differently, because they are different things. The first is a core, found from the gauge-invariant phase winding around a plaquette; the second has no core to find and is read from the winding on a contour drawn in the metal around the hole, which is an exact integer however little field threads the opening.

Flux expulsion by an S/I/S ring

A 1 µm hole centred in a 4 µm S/I/S plane with 500 nm layers, at ξ = 100 nm — Nb near $T_c$, where Ginzburg–Landau applies. The device expels flux completely — no vortices anywhere, zero fluxoid through the hole — up to 9.2 ± 0.3 mT.

What limits it is not the hole. A 4 µm plane is 20 λ across and screens so well that only 1.7% of the applied field reaches the hole (0.07 Φ₀ through it at the threshold), so the ring is nowhere near its fluxoid limit; vortices penetrate the 1.5 µm-wide arms first, and the hole does not admit a fluxoid until 10.9 mT. The device therefore beats the naive single-loop estimate $\Phi_0/A_\text{hole} = 2.07$ mT by more than a factor of four.

Nucleation and screening currents

Circuits

A superconducting bridge rectifier

The vortex diode of Castellani, Medeiros et al. (2024) is a wire with a triangular notch on one edge: in a perpendicular field the notch lets vortices in more easily for one current direction than the other. tdgl3d rebuilds the diode and then the paper’s full-wave rectifier — four diodes and a normal-metal load — as a single TDGL device.

The diode reproduces the measurement’s mechanism: equal critical currents at zero field (exact, by symmetry) and an efficiency rising linearly with field. The field scale is within a factor of 4–5 of the chip’s, and a static thin-film London solve of the wire between its contact pads accounts for that factor if the wire is short against pads of order 100 µm. In the bridge, nothing reaches the load below twice the reverse diodes’ critical current. Between that and twice the forward one, the load voltage has the same sign for either input polarity: under a sine drive the output is full-wave. Its average comes within 2% of the DC transfer curve averaged over the same sine.

Self-heating: a hotspot that switches and latches

Moving vortices and normal current dissipate energy. The dissipation density $Q = 2|(\partial_t + i\mu)\psi|^2 + 2|E|^2$ drives a heat equation for the reduced temperature θ, and θ lowers the local $T_c$ margin $\varepsilon = \varepsilon_0 - \theta$. In a strip with a notch, ramping the current 0 → 5 → 0, phase slips at the notch heat it past recovery at I ≈ 4.6, and a normal belt spreads across the width, sustained by its own Joule heat. On the way down it holds until I ≈ 1.9: the self-heating hysteresis of real nanowires. With heating off, the same ramp shows only a small phase-slip voltage and never switches.

Single-photon detection

A photon absorbed in a current-carrying strip heats a spot a few ξ across. The replication of Zotova and Vodolazov (2012) runs their NbN strip at T_c/2: the hot spot heals to |ψ| ≈ 0.3, a vortex–antivortex pair nucleates in it about 4 ps after absorption, the pair runs to the edges, and the heat it dissipates grows a normal belt that latches under the fixed bias. The one threshold the paper states — detection at 0.90 of the depairing current for a 52 ξ strip — comes out at 0.921, and 0.930 at half the grid spacing. The trends follow: the detection current falls with photon energy, edge absorption is easier than central at low energy and harder at high, and the threshold rises with width and saturates.

The same harness, scaled to a wide WSi strip after Patel et al. (JPL), measures latency against where the photon lands. From 0.25 to 4 µm wide at 0.8 I_dep, the time to a detectable voltage grows with width, and in strips 2 µm and wider a photon near the edge is detected sooner than one in the centre, as the measurement shows.

A detector block you can run

The live simulator runs the same 3-D scheme in the browser, on the reader’s GPU: a block of NbN on 20 nm of SiO₂ over silicon at 4.2 K, between two gold contacts, with ψ, the vector and scalar potentials, the current’s own field, and the heat equation in the film and the oxide. Drag to orbit it, set the bias, field and penetration depth, and click the film to absorb a photon. The film is coloured by the order parameter |ψ| or its phase, so vortices show as dark spots or as points where the phase colours wind. Traces show the current, voltage and peak temperature, and a ramp mode sweeps the bias from positive to negative to draw the full I–V curve, waiting at zero until the block is cold again. The page checks its GPU against a double-precision reference that matches tdgl3d to round-off before it runs.

The I–V shows why thick wires latch. The block switches near the depairing current straight onto its normal-state line, about 114 Ω across the probes for 39 nm, and then its own Joule heat holds it far above T_c. Heat leaves only through the oxide, and a 39 nm block takes about 280 ps to shed it against about 50 ps for a 5 nm film, so it stays normal until the current is nearly gone.

Vortices alone do not latch it. The Flux flow preset puts the 20 nm block in 0.5 H_c2 at a fifth of the depairing current: about ten vortices stream across and are replaced as they leave, the flux-flow voltage settles near a tenth of the normal value, and the film stays below about 5 K because the oxide carries that small, spread-out dissipation away. A little more bias or field and the vortices merge into a self-heated normal state, which is the runaway the ramp shows. A switch turns the Joule heating off for the isothermal flow of the 2-D demos.

Heater nanocryotrons

An hTron is a superconducting channel switched by heat from a resistive heater that is electrically isolated from it: the heater’s Joule heat crosses an insulator, lowers the channel’s switching current below its bias, and the bias current moves into the load. tdgl3d models the heater as a normal-metal region, the insulator as a dielectric that passes heat but no current, and puts the channel in a readout circuit with a bias, a load resistor and a series inductance.

The load decides whether it resets. With R_L = 10 instead of 3 more of the bias stays in the resistive channel (0.82 I_b), its own Joule heat keeps the hot region normal after the heater turns off, and the device latches — the electrothermal latching of nanowire switches, set by how the load shares the current.

Stacked in 3-D — a 1.5 ξ channel, 1 ξ of oxide and a heater line crossing over it — the heat crosses the oxide and spreads along the channel, and the channel’s switching current falls from 3.35 with no heater to 0.25 at a heater current of 6. Against the mean temperature under the heater, the first points follow 1 − θ more closely than the Ginzburg–Landau depairing law.

Inductance, impedance and capacitance

A circuit needs its lines’ L and C, and three solvers provide them:

• From TDGL — inductance per unit length is the stored field plus kinetic energy over $I^2$, and a small sinusoidal drive gives the complex impedance from the Poynting balance on one cross-section. In a field, the vortices add a Campbell-like inductance and dominate the loss at low frequency.
• Static London extraction — the port inductance matrix of 3-D superconducting wiring, with ports as cuts through the conductor, like a junction. It matches Swihart’s parallel-plate inductance to 0.39%, and two stacked 10 µm Nb rings come out at 19.9 pH each with k = 0.46.
• Electrostatics — a 2-D Laplace solve gives the capacitance matrix per unit length, matching closed forms for stripline and coplanar waveguide to 0.1–0.3% and Hammerstad and Jensen’s microstrip to 0.25%. With the TDGL inductance, that gives $Z_0 = \eta_0\sqrt{L/C}$ and the phase velocity, kinetic inductance included.

Checks against exact solutions

Two limits of the coupled equations have closed-form solutions, and between them they exercise each equation on its own. Both comparisons have no fitted parameters, and both run at three grid spacings so the residual can be shown to be discretisation error rather than disagreement.

In the London limit (|ψ| = 1, so the ψ-equation drops out) a square with the field pinned on its boundary obeys $\nabla^2 B = B/\lambda^2$, which has an exact Fourier solution. The solver matches it to rms 4.1e-3 · B₀ at h = 1 ξ, falling to 3.3e-4 at h = 0.25 ξ — observed order 1.82 in h. On the pinned boundary plaquettes, where the condition is Dirichlet and the solver should be exact rather than approximate, it agrees with the applied field to 6e-16.

At a pair-breaking wall (zero field, so the gauge field drops out) $\psi'' = -\psi + \psi^3$ gives $\tanh((x - x_0)/\sqrt{2})$, with the offset $x_0$ fixed by matching to the insulator’s relaxation rather than fitted. The solver matches it to rms 5.0e-2 at h = 1 ξ, falling to 4.8e-3 at h = 0.25 ξ — observed order 1.69 in h. The √2 is the physics being checked: the Ginzburg–Landau healing length is √2 ξ, not ξ.

The 3D path is checked against the same exact solution: a problem with no z-dependence must be solved identically by the 2D and 3D codes, and it is — the field varies across z-slices by 2e-16 and differs from the 2D run by 2e-10. Altogether the verification suites record 294 physics checks — gauge covariance, ∇·B = 0, symmetry, fluxoid quantisation, the closed-form limits — each with its measured value, the value physics requires, and the tolerance allowed.

Scale

The 3×3 hole array above is 1.8 M nodes at ξ = 100 nm, and it is the size that decides whether a study is an afternoon or a month. Measured on 4 cores, per unit of Ginzburg–Landau time, with forward Euler at 0.9 of the CFL limit:

All eight were run, not extrapolated. Three knobs matter at that scale: a thread pool for the right-hand side (bandwidth-bound, so cores help — 2.7× on four), streaming frames to HDF5 as they are produced so memory holds one frame however long the run is, and single precision, which cuts both the memory and the bandwidth the evaluation is limited by.

Where the field evolves, use the IMEX integrator. It takes the stiff $\kappa^2\nabla\times\nabla\times\mathbf{A}$ term implicitly, solved exactly by sine and cosine transforms, so the step is set by ψ rather than by κ: 11× faster than forward Euler at κ = 5 and 170× at κ = 20, agreeing with it to 1e-3 in ψ. Where A is frozen, forward Euler is the same thing. The implicit trapezoidal integrator is not worth it: its unpreconditioned Newton–GCR solve costs more than the larger step buys.

Getting started

git clone https://github.com/omedeiro/nanowire_tdgl.git
cd nanowire_tdgl/packages/tdgl3d
pip install -e ".[dev]"
pytest

import tdgl3d

params = tdgl3d.SimulationParameters(
    Nx=20, Ny=20, Nz=4,
    hx=1.0, hy=1.0, hz=1.0,
    kappa=5.0,
)
field = tdgl3d.AppliedField(Bz=1.0, ramp=True, ramp_fraction=0.3)
device = tdgl3d.Device(params, applied_field=field)

solution = tdgl3d.solve(device, t_stop=10.0, dt=0.05, method="euler")

solution.plot_order_parameter(slice_z=2)

The hole-array figures are reproduced by packages/tdgl3d/examples/nb_hole_array.py, whose --dry-run prints the grid, the memory per frame and a wall-time estimate before you commit to a run. Every figure above is produced by a standalone script in docs/figures/.

MATLAB predecessor

The Python package is a rewrite of the 3D TDGL MATLAB code written for MIT 6.336 (Spring 2021), verified against it index for index.

Source

github.com/omedeiro/nanowire_tdgl — the Python solver, the project schema and the job server. The original MATLAB is in github.com/omedeiro/simulation6336.