# SCRIPT: BH-BELL-11
# Bell at the horizon, two routes, each checked before it is written.
# Route E (entanglement): a Hawking pair of one frequency is a two-mode squeezed state
#   |psi> = sum_n c_n |n>_out |n>_in, c_n = tanh(r)^n / cosh(r), tanh(r) = exp(-omega / (2 T_H)).
# Route A (anomaly): eliminating the modes that fall in leaves a one-way outside theory whose balance law fails;
#   the Hawking flux is what restores it (Robinson and Wilczek 2005).
# Kill conditions, written before the run:
#  E1  the outside half alone is exactly thermal: its reduced state is diagonal with Bose weights tanh(r)^(2n)/cosh(r)^2,
#      and it is the same for every value of a correlation phase put on the pair (the local record cannot see it).
#  E2  with the parity pseudospin operators of Chen, Pan, Hou and Zhang, the CHSH value of the pair is
#      2 sqrt(1 + F^2) with F = tanh(2r) > 0: above the local-hidden-record bound 2 at every frequency. Kill if any
#      computed value differs from the formula by more than 1e-8 at the truncation used, or does not exceed 2.
#  E3  same local record, opposite futures (Atlas B108): (|10> +- i|01>)/sqrt2 have identical one-mode states and,
#      under the exchange g(|10><01| + |01><10|), opposite imbalances +-sin(2 g t).
#  A1  reproduction: with f(r) = 1 - 2M/r, surface gravity kappa = f'(r_H)/2 and T_H = kappa/(2 pi), the anomaly
#      cancellation flux f'(r_H)^2/(192 pi) equals the thermal flux pi T_H^2/12 of one chiral 1+1 boson. Symbolic.
#  A2  on our side the anomaly-shaped term is the re-entry of what a reduction discards: for x * y = P(xy),
#      (x*y)*z - x*(y*z) = P[x Q(yz) - Q(xy) z] (Atlas B105), and for the transport, P D^2 P - (P D P)^2 = P D Q D P.
#      Random checks to 1e-12.
#  A3  the one-way balance: in the retention model with return rate b = 0, the outside energy obeys
#      dx/dt = P - kappa x - a x exactly, so the outside books close only when the one-way flux a x is counted.
import numpy as np, sympy as sp
from scipy.linalg import expm
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
N=160
def pair_state(r,phase=0.0):
    t=np.tanh(r); c=np.array([t**n for n in range(N)])/np.cosh(r)
    c=c*np.exp(1j*phase*np.arange(N))
    psi=np.zeros((N,N),dtype=complex); psi[np.arange(N),np.arange(N)]=c
    return psi   # psi[m,n] amplitude of |m>_out |n>_in
# E1
r=0.7; rho_out=lambda p: p@p.conj().T
d0=rho_out(pair_state(r)); d1=rho_out(pair_state(r,1.3))
bose=np.array([np.tanh(r)**(2*n) for n in range(N)])/np.cosh(r)**2
offd=np.abs(d0-np.diag(np.diag(d0))).max()
report('E1 outside half is exactly thermal and blind to the pair phase', offd<1e-14 and np.abs(np.diag(d0).real-bose).max()<1e-14 and np.abs(d0-d1).max()<1e-14, f'off-diagonal {offd:.1e}, phase change {np.abs(d0-d1).max():.1e}')
# E2 pseudospin
sz=np.diag([(-1)**(n+1) for n in range(N)]).astype(complex)        # +1 on odd, -1 on even
sm=np.zeros((N,N),dtype=complex)
for k in range(N//2): sm[2*k,2*k+1]=1
sx=sm+sm.conj().T; sy=-1j*(sm-sm.conj().T) if False else 1j*(sm-sm.conj().T)
def corr(A,B,psi):
    # <psi| A (x) B |psi> with psi as matrix: sum conj(psi) * (A psi B^T)
    return np.sum(psi.conj()*(A@psi@B.T)).real
def chsh_max(psi):
    # maximise over spin directions: build the 3x3 correlation matrix T_ij and use Horodecki: 2 sqrt(sum of two largest squares)
    S=[sx,sy,sz]; T=np.array([[corr(a,b,psi) for b in S] for a in S])
    ev=np.sort(np.linalg.eigvalsh(T.T@T))[::-1]; return 2*np.sqrt(ev[0]+ev[1]), T
rows=[]
for r in [0.1,0.3,0.7,1.2]:
    val,T=chsh_max(pair_state(r)); F=np.tanh(2*r); pred=2*np.sqrt(1+F**2)
    rows.append((r,val,pred)); print(f'   r={r}: CHSH = {val:.10f}  formula 2 sqrt(1+tanh(2r)^2) = {pred:.10f}')
report('E2 CHSH above 2 at every squeezing, matching the formula', all(abs(v-p)<1e-8 and v>2 for _,v,p in rows))
# Hawking: tanh r = exp(-omega/2T): CHSH as a function of omega/T
for x in [0.5,1,2,5,10]:
    r=np.arctanh(np.exp(-x/2)); print(f'   omega/T_H = {x:>4}: squeezing r = {r:.4f}, CHSH = {2*np.sqrt(1+np.tanh(2*r)**2):.6f}')
# E3
g,t=0.7,0.4
def evolve(psi):
    H=g*(np.outer([0,1,0,0],[0,0,1,0])+np.outer([0,0,1,0],[0,1,0,0]))   # basis |00>,|10>,|01>,|11>
    return expm(-1j*H*t)@psi
p1=np.array([0,1,1j,0])/np.sqrt(2); p2=np.array([0,1,-1j,0])/np.sqrt(2)
def red(p):
    m=p.reshape(2,2)   # index [a][b] with |10> -> a=1,b=0 ; here ordering |ab>: |00>,|10>,|01>,|11> -> a = idx%2? use explicit
    return None
def local(p):
    # amplitudes: |00>=p0, |10>=p1 (mode1 excited), |01>=p2, |11>=p3
    M=np.array([[p[0],p[2]],[p[1],p[3]]])   # rows: mode-1 occupation, cols: mode-2 occupation
    return M@M.conj().T, M.T@M.conj()
l1=local(p1); l2=local(p2)
imb=lambda p: abs(p[1])**2-abs(p[2])**2
i1,i2=imb(evolve(p1)),imb(evolve(p2))
report('E3 same local record, opposite futures', np.allclose(l1[0],l2[0]) and np.allclose(l1[1],l2[1]) and abs(i1+i2)<1e-12 and abs(abs(i1)-np.sin(2*g*t))<1e-12, f'imbalances {i1:+.6f}, {i2:+.6f}; sin(2gt) = {np.sin(2*g*t):.6f}')
# A1
rr,M=sp.symbols('r M',positive=True); f=1-2*M/rr; rH=2*M
fp=sp.diff(f,rr).subs(rr,rH); kappa=fp/2; TH=kappa/(2*sp.pi)
flux_anom=fp**2/(192*sp.pi); flux_th=sp.pi*TH**2/12
report('A1 anomaly flux equals the thermal flux (reproduction)', sp.simplify(flux_anom-flux_th)==0, f'flux = {sp.simplify(flux_anom)}, T_H = {sp.simplify(TH)}')
# A2
rng=np.random.default_rng(11); err=0
for _ in range(20):
    n=6; X=rng.normal(size=(n,3)); P=X@np.linalg.solve(X.T@X,X.T); Q=np.eye(n)-P
    x,y,z=(rng.normal(size=(n,n)) for _ in range(3)); st=lambda a,b:P@(a@b)
    lhs=st(st(x,y),z)-st(x,st(y,z)); rhs=P@(x@Q@(y@z)-Q@(x@y)@z)
    D=rng.normal(size=(n,n)); D=D+D.T
    err=max(err,np.abs(lhs-rhs).max(),np.abs(P@D@D@P-(P@D@P)@(P@D@P)-P@D@Q@D@P).max())
report('A2 associator and compression defect identities', err<1e-12, f'max error {err:.1e}')
# A3
xs,ys,Pw,k,a=sp.symbols('x y P kappa a',positive=True)
xdot=Pw-(k+a)*xs+0*ys
report('A3 one-way books close only with the flux a x counted', sp.simplify(xdot-(Pw-k*xs)+a*xs)==0, 'dx/dt = P - kappa x - a x')
print('ALL PASS' if ok else 'SOME FAILED')