1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5# Reproducible finite-lattice test of the coherent-to-diffusive claim.
  6SEED = 1273
  7J = 1.0
  8N = 81
  9center = N // 2
 10
 11def rhs(rho, gamma):
 12    # H=-J(sum |x><x+1| + h.c.), and pure dephasing of off-diagonal entries.
 13    H = np.zeros((N, N), dtype=np.complex128)
 14    q = np.arange(N - 1)
 15    H[q, q + 1] = H[q + 1, q] = -J
 16    return -1j * (H @ rho - rho @ H) - gamma * (rho - np.diag(np.diag(rho)))
 17
 18def run(gamma, T=10.0, dt=0.01):
 19    rho = np.zeros((N, N), dtype=np.complex128)
 20    rho[center, center] = 1.0
 21    times, msd, corr, hf = [], [], [], []
 22    steps = int(round(T / dt))
 23    for s in range(steps + 1):
 24        if s % 5 == 0:
 25            p = np.real(np.diag(rho)); p = np.maximum(p, 0)
 26            x = np.arange(N) - center
 27            msd.append(float(np.sum(x*x*p) / max(np.sum(p), 1e-12)))
 28            diag_norm = np.linalg.norm(np.real(np.diag(rho)))
 29            corr.append(float(np.linalg.norm(rho - np.diag(np.diag(rho))) / (diag_norm + 1e-12)))
 30            # high spatial frequencies of the population (proxy for oscillations)
 31            fft = np.abs(np.fft.rfft(p - p.mean())) ** 2
 32            hf.append(float(fft[len(fft)//2:].sum() / (fft.sum() + 1e-12)))
 33            times.append(s * dt)
 34        if s == steps: break
 35        # RK4 keeps the integration explicit and makes the mechanism auditable.
 36        k1 = rhs(rho, gamma)
 37        k2 = rhs(rho + .5*dt*k1, gamma)
 38        k3 = rhs(rho + .5*dt*k2, gamma)
 39        k4 = rhs(rho + dt*k3, gamma)
 40        rho = rho + dt*(k1 + 2*k2 + 2*k3 + k4)/6
 41        rho = (rho + rho.conj().T) / 2
 42        rho /= np.trace(rho).real
 43    return np.array(times), np.array(msd), np.array(corr), np.array(hf)
 44
 45def fit_power(t, y, lo=3.0):
 46    z = (t >= lo) & (y > 1e-10)
 47    return float(np.polyfit(np.log(t[z]), np.log(y[z]), 1)[0])
 48
 49def fit_slope(t, y, lo=5.0):
 50    z = t >= lo
 51    return float(np.polyfit(t[z], y[z], 1)[0])
 52
 53def feedback(T=10.0, dt=.01, gamma0=.5, target=.10, alpha=.8, gmin=.02, gmax=32.):
 54    rho = np.zeros((N,N), complex); rho[center,center] = 1
 55    x=np.arange(N)-center; t=[]; m=[]; gs=[]; rs=[]
 56    gamma=gamma0; steps=int(T/dt)
 57    for s in range(steps+1):
 58        if s % 5 == 0:
 59            p=np.maximum(np.real(np.diag(rho)),0)
 60            r=np.linalg.norm(rho-np.diag(np.diag(rho)))/(np.linalg.norm(np.real(np.diag(rho)))+1e-12)
 61            t.append(s*dt); m.append(np.sum(x*x*p)/max(p.sum(),1e-12)); gs.append(gamma); rs.append(r)
 62        if s==steps: break
 63        k1=rhs(rho,gamma); k2=rhs(rho+.5*dt*k1,gamma); k3=rhs(rho+.5*dt*k2,gamma); k4=rhs(rho+dt*k3,gamma)
 64        rho=(rho+dt*(k1+2*k2+2*k3+k4)/6); rho=(rho+rho.conj().T)/2; rho/=np.trace(rho).real
 65        # controller is stopped-gradient by construction: r is a scalar observation.
 66        if s % 5 == 4:
 67            gamma=float(np.clip(gamma*np.exp(alpha*(r-target)),gmin,gmax))
 68    return np.array(t),np.array(m),np.array(gs),np.array(rs)
 69
 70def main():
 71    np.random.seed(SEED)
 72    gammas=[0., .5, 1., 2., 4., 8., 16.]
 73    rows=[]
 74    for g in gammas:
 75        t,m,r,h=run(g)
 76        slope=fit_slope(t,m); beta=fit_power(t,m)
 77        # D prediction says MSD slope = 2D = 4J^2/gamma.
 78        pred=None if g==0 else 4*J*J/g
 79        rows.append(dict(gamma=g, msd_final=float(m[-1]), late_msd_slope=slope,
 80                         predicted_slope=pred, slope_ratio=None if pred is None else slope/pred,
 81                         power_exponent=beta, corr_final=float(r[-1]), hf_final=float(h[-1])))
 82    # Selected-frequency crossover: D k^2 = J |sin k|, gamma=2J k^2/|sin k|.
 83    kstar=np.pi/2; predicted_cross=2*J*kstar*kstar/abs(np.sin(kstar))
 84    # empirical crossover: nearest beta=1.5, with gamma~4J prediction separately.
 85    empirical=min(rows[1:], key=lambda z: abs(z['power_exponent']-1.5))['gamma']
 86    tc,mm,gg,rr=feedback()
 87    out={
 88      'parameters': {'J':J,'N':N,'T':10.,'dt':.01,'seed':SEED},
 89      'predictions': {
 90        'diffusion_slope': 'MSD/t = 4 J^2/gamma',
 91        'crossover_bandwidth': 'gamma approximately 4J',
 92        'selected_k_crossover': predicted_cross,
 93        'ballistic_exponent': 2.0, 'diffusive_exponent': 1.0
 94      },
 95      'sweep': rows,
 96      'observed': {'beta_1p5_crossover_gamma': empirical,
 97                   'feedback_initial_gamma':.5,'feedback_final_gamma':float(gg[-1]),
 98                   'feedback_final_msd':float(mm[-1]),'feedback_max_gamma':float(gg.max())}
 99    }
100    Path('results.json').write_text(json.dumps(out, indent=2))
101    print(json.dumps(out, indent=2))
102
103if __name__=='__main__': main()