# SCRIPT: BH-OFFSET-15
# The offset of the outside half of a fermionic Hawking pair.
# State of one pair (Dirac mode of frequency w, Hawking temperature T): cos r|00> + e^{i phi} sin r|11>,
# tan^2 r = exp(-w/T) (Fermi-Dirac). Outside covariance C = <c^dagger c> on the outside mode(s).
# Offset (Beasley, The Offset Belongs to the Boundary, v31): tau = log det(I-C)/det C = log P(empty)/P(full).
# Why this tests the claim: the outside reduced state is computed from the full two-mode state vector by a
# partial trace (no thermal formula is assumed), so if the outside record were not exactly Fermi-Dirac,
# or if the pair phase leaked into it, tau would differ from w/T or depend on phi.
# Kill conditions, written before the run:
#  K1  one pair: the outside reduced density matrix is diagonal for every phi (off-diagonal < 1e-15), and
#      tau = w/T to 1e-13 over six frequencies and seven phases. Kill otherwise.
#  K2  four independent pairs (product state, 4 outside modes, 16-dimensional outside Fock space): tau computed
#      from the 4x4 outside covariance equals log P(all empty)/P(all full) from the Fock-space reduced state,
#      and both equal sum_k w_k/T, to 1e-12, with random phases. Kill otherwise.
#  K3  scale: replacing (w, T) by (s w, s T) leaves tau unchanged (it is the ratio w/T), so the offset does not
#      supply an absolute scale. Kill if it changes by more than 1e-13.
import numpy as np, itertools
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
def pair_state(w,T,phi):
    r=np.arctan(np.exp(-w/(2*T)))  # tan r = exp(-w/2T), so tan^2 r = exp(-w/T)
    psi=np.zeros(4,complex); psi[0]=np.cos(r); psi[3]=np.exp(1j*phi)*np.sin(r)  # basis |out,in> = 00,01,10,11
    return psi
def outside_rho(psi):
    M=psi.reshape(2,2); return M@M.conj().T
err=0; offd=0
for w in [0.1,0.5,1.0,2.0,5.0,12.0]:
    for phi in np.linspace(0,2*np.pi,7,endpoint=False):
        rho=outside_rho(pair_state(w,1.0,phi)); offd=max(offd,abs(rho[0,1]))
        C=rho[1,1].real; tau=np.log((1-C)/C); err=max(err,abs(tau-w))
report('K1 one pair: outside record Fermi-Dirac, tau = w/T for every phase', offd<1e-15 and err<1e-13, f'max off-diagonal {offd:.1e}, max |tau - w/T| {err:.1e}')
rng=np.random.default_rng(15); ws=np.array([0.7,1.3,2.9,8.36]); T=1.0
psi=np.array([1.0+0j])
for w in ws: psi=np.kron(psi,pair_state(w,T,rng.uniform(0,2*np.pi)))
# reorder to (out1..out4, in1..in4)
psi=psi.reshape([2]*8); psi=np.transpose(psi,[0,2,4,6,1,3,5,7]).reshape(16,16)
rho=psi@psi.conj().T
def occ(i,k): return (i>>(3-k))&1
Cmat=np.zeros((4,4),complex)
for i,j in itertools.product(range(16),range(16)):
    for a in range(4):
        for b in range(4):
            # <c_a^dagger c_b> for a != b needs a Fock off-diagonal element of rho; those are checked to vanish below (fock_offd)
            if a==b and i==j and occ(i,a): Cmat[a,b]+=rho[i,i]
offC=np.abs(Cmat-np.diag(np.diag(Cmat))).max()
tauC=np.log(np.linalg.det(np.eye(4)-Cmat).real/np.linalg.det(Cmat).real)
tauF=np.log(rho[0,0].real/rho[15,15].real)
fock_offd=np.abs(rho-np.diag(np.diag(rho))).max()
report('K2 four pairs: tau from covariance = log P(empty)/P(full) = sum w_k/T', abs(tauC-ws.sum())<1e-12 and abs(tauF-ws.sum())<1e-12 and fock_offd<1e-15,
       f'tau_C = {tauC:.12f}, tau_Fock = {tauF:.12f}, sum w/T = {ws.sum():.12f}; outside Fock off-diagonal {fock_offd:.1e}')
d=0
for s in [0.01,3.0,1e3]:
    rho=outside_rho(pair_state(s*2.0,s*1.0,0.4)); C=rho[1,1].real; d=max(d,abs(np.log((1-C)/C)-2.0))
report('K3 the offset is the ratio w/T; a common rescaling leaves it unchanged', d<1e-13, f'max change {d:.1e}')
print('ALL PASS' if ok else 'SOME FAILED')