Yukawa monolayer -- 2D dusty plasma MD

—

integration unstable — reduce the time step

Run

Confinement

In-plane potential of a ring of radius a lying h below the layer. Matched to the parabola at the center, so it differs only where the cluster reaches the ring.

Ring radius a 4.50

Ring depth η = h/a 0.25

Confines only for η < 0.707. Deeper rings give a shallower, softer well.

Plasma conditions

Grains N 100

Changing N reseeds the run.

Screening κ = r₀/λ_D 2.00

0 is bare Coulomb. Higher κ shortens the range and shrinks the cluster.

Temperature T* 1.0e-4

Langevin setpoint, in units of mω₀²r₀²/k_B.

Gas drag ν 0.50

Epstein friction. At ν = 0 the thermostat switches off and the run is microcanonical — watch the energy trace stay flat.

Display

Color grains by

Delaunay bonds Radial grid (1 r₀) Motion trails

Numerics

Time step Δt 0.020

Steps per frame 8

Diagnostics

Model, units and verification

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.
ψ₆ 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 ψ₆ = 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.