1import json
 2import numpy as np
 3
 4SX=np.array([[0,1],[1,0]],complex)
 5SY=np.array([[0,-1j],[1j,0]],complex)
 6SZ=np.array([[1,0],[0,-1]],complex)
 7I2=np.eye(2,dtype=complex)
 8
 9def H(kx,ky,m,q=1.):
10    return q*np.sin(kx)*SX + q*np.sin(ky)*SY + (m+np.cos(kx)+np.cos(ky))*SZ
11
12def drift(kx,ky,m,gamma,q=1.):
13    return -gamma*I2-1j*H(kx,ky,m,q)
14
15def spectrum(kx,ky,m,gamma,omega=0.8,noise=1.,q=1.):
16    A=drift(kx,ky,m,gamma,q)
17    G=np.linalg.inv(-1j*omega*I2-A)
18    return G@(noise*I2)@G.conj().T
19
20def link(a,b):
21    z=np.vdot(a,b)
22    return z/(abs(z)+1e-30)
23
24def fukui_chern(m,gamma=.7,omega=.8,n=41,band=0):
25    us=np.empty((n,n,2),complex)
26    for ix in range(n):
27      for iy in range(n):
28        w,v=np.linalg.eigh(spectrum(2*np.pi*ix/n,2*np.pi*iy/n,m,gamma,omega))
29        us[ix,iy]=v[:,::-1][:,band]
30    total=0.
31    for ix in range(n):
32      for iy in range(n):
33        u=us[ix,iy]; ux=us[(ix+1)%n,iy]; uy=us[ix,(iy+1)%n]
34        uxy=us[(ix+1)%n,(iy+1)%n]
35        plaquette=link(u,ux)*link(ux,uxy)*link(uxy,uy)*link(uy,u)
36        total += np.angle(plaquette)
37    return total/(2*np.pi)
38
39def min_gap(m,gamma=.7,omega=.8,n=81,noise=1.):
40    gaps=[]
41    for kx in np.linspace(-np.pi,np.pi,n,endpoint=False):
42      for ky in np.linspace(-np.pi,np.pi,n,endpoint=False):
43        w=np.linalg.eigvalsh(spectrum(kx,ky,m,gamma,omega,noise))
44        gaps.append(w[1]-w[0])
45    return float(np.min(gaps))
46
47def max_h2(m,n=81):
48    return max(np.linalg.eigvalsh(H(x,y,m)).max()**2 for x in np.linspace(-np.pi,np.pi,n,endpoint=False) for y in np.linspace(-np.pi,np.pi,n,endpoint=False))
49
50def euler_radius(dt,m,gamma,n=41):
51    r=0.
52    for x in np.linspace(-np.pi,np.pi,n,endpoint=False):
53      for y in np.linspace(-np.pi,np.pi,n,endpoint=False):
54        ev=np.linalg.eigvals(I2+dt*drift(x,y,m,gamma))
55        r=max(r,float(np.max(np.abs(ev))))
56    return r
57
58def strip_edge(m,gamma=.18,nx=40,ny=20,q=1.):
59    best=[]
60    chi=q
61    for kx in 2*np.pi*np.arange(nx)/nx:
62      M=np.zeros((2*ny,2*ny),complex)
63      for y in range(ny):
64        M[2*y:2*y+2,2*y:2*y+2]=(m+np.cos(kx))*SZ+chi*np.sin(kx)*SX
65        if y+1<ny:
66          T=(-1j*chi*SY+SZ)/2
67          M[2*y:2*y+2,2*(y+1):2*(y+1)+2]=T
68          M[2*(y+1):2*(y+1)+2,2*y:2*y+2]=T.conj().T
69      w,v=np.linalg.eigh(M)
70      for j,e in enumerate(w):
71        if abs(e)<0.35:
72          p=np.sum(abs(v[:,j].reshape(ny,2))**2,axis=1)
73          edge=(p[0]+p[1]+p[-1]+p[-2])/(p.sum()+1e-15)
74          ipr=float(np.sum(p*p))
75          best.append((edge,ipr,float(e)))
76    return max(best) if best else (0.,0.,999.)
77
78def main():
79  np.random.seed(0); gamma=.7; omega=.8
80  lam=max_h2(-1.); dtcrit=2*gamma/(gamma*gamma+lam)
81  dts=[.8*dtcrit,.99*dtcrit,1.01*dtcrit,1.2*dtcrit]
82  stability=[(dt,euler_radius(dt,-1.,gamma)) for dt in dts]
83  masses=[-2.5,-1.5,-.5,.5,1.5,2.5]
84  phase=[(m,round(fukui_chern(m,gamma,omega,n=31)),min_gap(m,gamma,omega,n=51)) for m in masses]
85  near=[(m,round(fukui_chern(m,gamma,omega,n=31),3),min_gap(m,gamma,omega,n=61)) for m in [-.2,-.1,-.05,.05,.1,.2]]
86  edge_top=strip_edge(-1.,q=1.); edge_triv=strip_edge(-1.,q=0.)
87  chirality=[(q,strip_edge(-1.,nx=32,ny=16,q=q)[:2]) for q in [0.,.1,.25,.5,1.]]
88  exact_gap=[(m,min_gap(m,gamma,omega,n=81)) for m in [-2.,0.,2.]]
89  out={'predictions':{
90    'stability_bound':{'predicted_dt_critical':dtcrit,'observed_dt_radius_1_crossing':stability,'pass': bool(all((r<1)==(dt<dtcrit) for dt,r in stability))},
91    'chern_gap_transition':{'predicted_boundaries':[-2,0,2],'phase_sweep':phase,'near_m0':near,'exact_boundary_gaps':exact_gap,'pass': bool(abs(phase[1][1])==1 and abs(phase[-2][1])==1 and min(x[2] for x in near)<max(x[2] for x in phase))},
92    'boundary_localization':{'topological_m_minus1_edge_ratio_ipr':edge_top,'deterministic_nonchiral_q0_baseline':edge_triv,'chirality_sweep':chirality,'pass':edge_top[0]>edge_triv[0]*2 and edge_top[1]>edge_triv[1]*1.3}
93  },'parameters':{'gamma':gamma,'omega':omega,'grid':31}}
94  print(json.dumps(out,indent=2))
95if __name__=='__main__': main()