What is being solved
N grains of charge Q and mass m, confined to a plane by an isotropic parabolic
trap and repelling through a Debye-screened Coulomb potential:
U_ij = Q²/(4πε₀r) · exp(-r/λ_D)
U_trap = ½ m ω₀² r²
Neutral gas supplies Epstein drag at rate ν and, through the
fluctuation–dissipation theorem, the matching Langevin noise.
Reduced units
length r₀ = [2Q²/(4πε₀mω₀²)]^(1/3)
time 1/ω₀
energy mω₀²r₀²
r₀ is exactly the equilibrium separation of two unscreened grains in
the trap, which gives the code a parameter-free analytic benchmark. The equation of
motion becomes
d²x/dt² = -x
+ ½ Σ (1/r² + κ/r) e^(-κr) rhat
- νv + noise
Confinement
The layer is strictly two-dimensional — there is no vertical trap. What confines the
grains is the in-plane potential of a metallic ring of radius a lying
a distance h below the layer, which raises the potential energy along that
annulus. Averaged over azimuth,
W(r) = 2/(π√((a+r)²+h²)) · K(m)
m = 4ar/((a+r)²+h²)
with K the complete elliptic integral of the first kind, evaluated by AGM. The amplitude
is fixed by U''(0) = 1, so the ring and the parabola share their central
curvature and r₀ keeps its meaning in both.
Two things differ from a parabola: the well has finite depth, and near the crest
at r ≈ a the restoring force is weaker, so the outer shell pushes further
out. Both only matter once the cluster grows out to meet the ring. Confinement requires
η = h/a < 1/√2; beyond the crest the barrier is capped so grains cannot escape
the box, while inside it the potential is exact.
Integrator
Grønbech-Jensen–Farago Langevin–Verlet (GJF-2GJ). At ν = 0
it reduces exactly to velocity Verlet, so turning the drag off gives a genuine
microcanonical run you can use as an energy-conservation test.
Temperature is read from the 2GJ half-step velocity, which samples the kinetic
energy without time-step bias.
Diagnostics
- Γ = 1/(2 b T*), with b the mean nearest-neighbor distance measured
from the Delaunay graph. This is the bare Coulomb coupling; the screened value
Γe-κb is also reported.
- ψ6 is the local hexagonal bond-order parameter over Delaunay
neighbors. It runs to 1 in a perfect triangular lattice.
- Defects are interior grains whose coordination is not 6. Boundary grains
are excluded, since they are undercoordinated by construction.
Verified against
- N = 2, κ = 0 → separation 1.000000 (exact, defines r₀)
- N = 3 → ring radius 0.660901 vs analytic 0.660901
- N = 4 → ring radius 0.782186 vs analytic 0.782186
- NVE drift |δE/E| < 10-8 over 50k steps at Δt = 0.02
- Thermostat reproduces the setpoint to within 1%
- Delaunay on a perfect lattice: Euler relation exact, interior ψ6 = 1.000000
Known limits
- Forces are all-pairs, O(N²). Fine to a few hundred grains; a neighbor list
is needed beyond that.
- Strictly 2D. No out-of-plane motion, so the ion-wake-driven vertical
instability that drives real monolayer melting is absent.
- The trap is a fixed parabola; grain charge does not vary with position.
- Shell occupations are read off a time-averaged profile and will flicker near
a structural transition.