Build a single self-contained HTML file (no build step, external assets limited to a plotting library and/or the Three.js CDN if needed) that implements and visualizes a Monte Carlo simulation of compact U(1) lattice gauge theory in 2+1 dimensions (pure gauge, no matter fields). The physics must be correct and checkable against known results. Match the following specifications.
1. Lattice & Field Definitions

A 3D periodic cubic lattice of size L³ (default L=8, adjustable 4–16). Dimensions are Euclidean (2 space + 1 Euclidean time); periodic boundary conditions in all three directions.
Degrees of freedom are compact U(1) link variables U_μ(n) = exp(iθ_μ(n)), with one angle θ_μ(n) ∈ (−π, π] per link, for each site n and direction μ ∈ {0,1,2}. Store the angles, not the complex numbers, to save memory.
The plaquette is the counterclockwise product of four links around an elementary square: U_μν(n) = U_μ(n) U_ν(n+μ̂) U_μ(n+ν̂)* U_ν(n)*, whose phase is θ_P = θ_μ(n) + θ_ν(n+μ̂) − θ_μ(n+ν̂) − θ_ν(n).

2. Action & Update Algorithm

Use the Wilson action S = β · Σ_P [1 − cos θ_P], summed over all distinct plaquettes. β = 1/g² is the inverse coupling (adjustable, default sweep range ~0.5 to 2.5).
Implement a correct local Metropolis update: for each link, propose θ → θ + δ with δ uniform in [−ε, ε], accept with probability min(1, exp(−ΔS)) where ΔS depends only on the (4 in 2+1D, i.e. the staples touching that link) plaquettes containing the link. Auto-tune ε to keep acceptance near ~50%, or expose ε as a control.
Optionally also offer a heat-bath update for U(1) as a faster alternative (stretch; Metropolis is the required baseline).
One "sweep" updates every link once. Support thermalization sweeps (discarded) followed by measurement sweeps with a configurable interval between measurements to reduce autocorrelation.

3. Required Observables (with known-answer checks)

Average plaquette ⟨cos θ_P⟩ as a function of β. This is the primary validation observable — at small β it should approach the strong-coupling result ⟨cos θ_P⟩ ≈ β/2, and at large β approach 1 − 1/(2·(something β-dependent)); display the measured curve and overlay the small-β analytic line so correctness is visually obvious.
Wilson loops W(R,T): rectangular R×T loops measured as the expectation of cos of the summed link phases around the loop. Measure for a range of R and T.
Static potential V(R), extracted from V(R) ≈ −(1/T) ln W(R,T) at the largest available T (and/or from ratios W(R,T)/W(R,T+1)). Plot V(R) vs R.
String tension via Creutz ratios: χ(R,T) = −ln[ W(R+1,T+1)·W(R,T) / ( W(R+1,T)·W(R,T+1) ) ]. Plot χ vs β. In 2+1D compact U(1) the theory is confining for all β (Polyakov), so a nonzero string tension that decreases with β is the expected signature — display this.

4. Visualization (this is the frontend test)

Live lattice view: a 2D slice (fix the third coordinate) showing plaquette energies as a heatmap — color each plaquette by its action density (1 − cos θ_P), e.g. cool color #1A3A6B for low energy → hot #F03E3E for high energy. The slice updates live as sweeps run, so the viewer sees fluctuations thermalize.
Live plots (use a CDN charting lib): (a) average plaquette vs sweep number (Monte Carlo history, should equilibrate to a plateau), (b) ⟨cos θ_P⟩ vs β scan with the analytic strong-coupling overlay, (c) static potential V(R) vs R, (d) Creutz-ratio string tension vs β.
A Wilson-loop diagram: let the user pick R and T with sliders and highlight that rectangular loop on the lattice slice, displaying its measured ⟨W(R,T)⟩ and the implied potential.

5. Controls (panel, live-updating)

Sliders/inputs for: lattice size L, β (single value), Metropolis step size ε, thermalization sweeps, measurement sweeps, measurement interval.
A "Run β scan" mode that sweeps β across a range and plots ⟨cos θ_P⟩ and the Creutz ratio across the scan (this is the headline confinement result).
Start / pause / reset. A live readout of: current sweep, acceptance rate, current ⟨cos θ_P⟩, and wall-clock sweeps/sec.

6. Performance & Correctness Constraints

The MC inner loop must be efficient (typed arrays for the angle field; compute ΔS from only the affected staples, never re-sum the whole action per link). Run the simulation in chunks via requestAnimationFrame (or a Web Worker — stretch) so the UI stays responsive and the lattice view animates.
The acceptance rate, the strong-coupling plaquette limit, and the equilibration plateau are the three things a correct implementation must reproduce; the prompt is considered failed if ⟨cos θ_P⟩ does not match β/2 at small β.
Single .html file. All JS inline. No localStorage/sessionStorage. No console errors. Code commented at each major section, especially the action, the staple/ΔS computation, and each observable.

7. Optional stretch modules (clearly separated, only if the core is solid)

SU(2) pure Yang-Mills: replace U(1) angles with SU(2) link matrices (a₀·𝟙 + ia·σ, |a|=1, stored as 4-vectors), Wilson action S = β·Σ_P[1 − ½Tr U_P], and a correct SU(2) heat-bath update drawing a₀ ∝ √(1−a₀²)·exp(βk·a₀) from the staple V (k = √det V). Validate against the known SU(2) average-plaquette curve. Note this is much heavier; keep the lattice tiny.
Flux-tube / confinement toy: visualize the action-density profile of the field between two static charges (Polyakov lines) to show the flux tube forming a linear-potential string.
Asymptotic freedom / running coupling: for the SU(N) case, plot the measured string tension or coupling vs β and compare to the perturbative two-loop running of g(β); show the coupling weakening at large β.

Success criteria for comparison: physical correctness above all — does ⟨cos θ_P⟩ reproduce the strong-coupling β/2 limit and equilibrate to a stable plateau; is the ΔS/staple computation correct (not a full-action recompute); is the plaquette/Wilson-loop geometry right (counterclockwise products, proper periodic wrapping); are the Creutz ratios and confinement signature genuinely emergent rather than hardcoded; quality and responsiveness of the live lattice heatmap and plots; and clean, well-commented code. Secondary: how far into the SU(2) and running-coupling stretch modules the model gets while staying correct.