Unexecuted reproducibility supplement

def synthetic_audit(): K = np.array([[.15,.05j],[-.05j,.15]],complex) rho,basis,tail,p,Cn = gaussian_cutoff(K,1,"synthetic_N1") ix = [basis.index((1,0)),basis.index((0,1))] expected = np.array([[425,125j],[-125j,425]],complex)/4356 check("synthetic_one_particle_entries",rho[np.ix_(ix,ix)],expected) check("synthetic_moment_orientation",Cn,expected.T) need(abs(Cn[0,1]-expected[0,1])>1e-3, "phase-sensitive regression is not discriminating") check("synthetic_eigenvalues",np.linalg.eigvalsh(K),np.array([.1,.2])) check("synthetic_vacuum",np.array(p[0]),np.array(25/33)) check("synthetic_tail",np.array(tail),np.array(103/2178)) gaussian_cutoff(K,3,"synthetic_N3") DATA["synthetic"] = "field_dataset":False,"a":8/33, "orientation_test":"explicit imaginary entries and truncated moments" dump() return 8/33

def zcoord(x): # Pilot |x| <= 80; retain z rather than forming r-1 afterwards. return float(lambertw(np.exp(x-1)).real)

def potential(x,ell): z = zcoord(x) r = 1+z return (z/r)*(ell*(ell+1)/r**2+1/r**3)

def radial_tail(cut,ell,omega_min): za,zb = zcoord(-cut),zcoord(cut) L = ell*(ell+1) v = L+.5 left = L*za/(1+za)+.5*za*(2+za)/(1+za)**2 right = L/(1+zb)+.5/(1+zb)**2 d = left+right w = v-d need(d>=0 and w>=0,"potential-tail quantities") A = d/(2*omega_min)*np.exp((v+w)/omega_min)*(1+w/omega_min) return "x_a":-cut,"x_b":cut,"r_h":1, "omega_min":omega_min,"left_tail":left,"right_tail":right, "d":d,"v":v,"w":w,"A":float(A), "S_error_bound_formula_value":float(min(2,4*A)), "vacuous_for_unitary_difference":bool(4*A>=2), "evaluation":"floating-point evaluation of analytical formula", "includes_ODE_error":False,"interval_certified":False

def radial(omega,ell,cut,rtol): def rhs(x,y): return np.array([y[1],(potential(x,ell)-omega**2)*y[0]],complex) y0 = np.array([np.exp(1j*omega*cut), 1j*omega*np.exp(1j*omega*cut)],complex) sol = solve_ivp(rhs,(cut,-cut),y0,method="DOP853", rtol=rtol,atol=rtol*1e-2) need(sol.success,"radial integration") f,fp = sol.y[:,-1] xa = -cut a = .5*np.exp(-1j*omega*xa)*(f+fp/(1j*omega)) b = .5*np.exp(1j*omega*xa)*(f-fp/(1j*omega)) tau,rh = 1/a,b/a ri = -np.conj(rh)*tau/np.conj(tau) S = np.array([[rh,tau],[tau,ri]],complex) flux = abs(abs(a)**2-abs(b)**2-1) unitary = np.linalg.norm(S.conj().T@S-I,2) return tau,rh,S,float(flux),float(unitary),int(sol.nfev)

def simpson_weights(grid): need(len(grid)>=3 and len(grid) weights = np.ones(len(grid)) weights[1:-1:2] = 4 weights[2:-1:2] = 2 return weights*(grid[-1]-grid[0])/(3*(len(grid)-1))