feat: add jump operator 𝒥f(x) with rigorous tests and README explainer

[Copilot]

Implements the Lévy–Khintchine nonlocal generator component — the only term in ℒ = b∂ₓ + ½σ²∂ₓ² + 𝒥 that evaluates f at the destination rather than derivatives at the origin.

New method: jump_operator

Monte Carlo evaluation of 𝒥f(x) = λ(x) ∫ [f(x+z) − f(x)] ν(dz|x):

eng = JumpDiffusionEngine(lambda_func, jump_rate=0.5, seed=42)

# Quadratic identity: 𝒥(x²) = λ for N(0,1) kernel, any x
eng.jump_operator(lambda x: x**2, x=2.0)  # → ~0.5

# Flat obstruction: all derivatives of e^{-1/x} vanish at 0, but 𝒥f(0⁺) > 0
eng.jump_operator(lambda x: np.exp(-1/x) if x > 1e-10 else 0.0, x=0.0)  # → > 0

Supports position-dependent jump_rate (callable) and custom jump_size_dist kernels. Deterministic by default via seed parameter.

Tests (TestJumpOperator, 9 cases)

• Exact: annihilates constants; zero when jump_rate=0
• Algebraic: linearity 𝒥(αf+βg) = α𝒥f+β𝒥g; quadratic identity 𝒥(x²)=λ; linear-function zero under symmetric kernel; rate scaling
• Structural: flat obstruction non-zero at x=0; crossing-condition cliff (P<0.15 below J_c, P>0.50 above, gap>0.35); custom kernel respected

The crossing condition test requires a bistable regime (k=0.5, g=2.0, K=1.0) to get well-defined basin walls.

README

Compact ### The Jump Operator 𝒥 subsection added directly under the Core Equation — formula, parameter mapping, confirmed properties, flat obstruction, and crossing condition in ~8 lines.

[beanapologist]

Add another example script

import numpy as np
from math import comb

print("FOCK SPACE — the occupation-number formulation")
print("="*70)

N=12 # truncated boson mode
a=np.diag(np.sqrt(np.arange(1,N)),1) # annihilation
ad=a.conj().T # creation
n_op=ad@a
I=np.eye(N)

print("1. THE LADDER OPERATORS (this is the jump operator, discretized)")
print(" [a,a†] =", np.round(np.diag(a@ad-ad@a)[:5],3), "... = 1 (canonical)")
print(" n̂|n> = n|n>, spectrum:", np.round(np.diag(n_op)[:8],1))
print(" a† JUMPS n->n+1. There is NO operator taking n->n+0.5.")
print(" -> Fock space IS the statement that transitions are jumps. Δn ∈ ℤ.")

print("\n2. SELECTION RULES = the jump kernel ν(dz|x), discretized")
print(" ⟨m|a†|n⟩ = √(n+1) δ_{m,n+1} -> kernel supported ONLY on z=+1")
print(" ⟨m|a |n⟩ = √n δ_{m,n-1} -> kernel supported ONLY on z=-1")
mat=np.round(ad[:5,:5],3)
print(" a† matrix (first 5x5):"); print(mat)
print(" every entry off the ±1 diagonal is ZERO: gradual transitions FORBIDDEN.")

print("\n3. COHERENT STATE = the classical limit (where jumps blur into a flow)")
for alpha in (0.5, 2.0, 6.0):
n_bar=alpha**2
dn=np.sqrt(n_bar) # Poisson: Δn=√n̄
print(f" |α|={alpha:>4}: n̄={n_bar:>6.1f} Δn={dn:>5.2f} Δn/n̄={dn/n_bar:>6.3f}"
+ (" <- discrete, jumps visible" if dn/n_bar>0.3 else " <- continuum limit"))
print(" Δn/n̄ = 1/√n̄ -> 0: THIS is why the wobble (n10^52) looks continuous")
print(" and J-space (n25) does not. Same operator, different n.")

print("\n4. TIME CRYSTAL IN FOCK LANGUAGE")
print(" DTC period-2: the Floquet operator U satisfies U²|ψ⟩=|ψ⟩ but U|ψ⟩≠|ψ⟩")

two-level cat: |ψ±> = (|0>+/-|1>)/sqrt2 swapped by U each period

psi_p=np.zeros(N); psi_p[0]=psi_p[1]=1/np.sqrt(2)
psi_m=np.zeros(N); psi_m[0]=1/np.sqrt(2); psi_m[1]=-1/np.sqrt(2)
U=np.eye(N); U[0,0]=1; U[1,1]=-1 # phase flip: swaps psi_p <-> psi_m
print(" U|ψ+⟩ = |ψ−⟩ , U|ψ−⟩ = |ψ+⟩ -> U²=1, period doubling")
print(" check U²=I:", np.allclose(U@U,np.eye(N)))
print(" the ORDER PARAMETER is which branch you're on: a BINARY, not an angle.")
print(" Melting = a jump between branches, not a rotation. Rigidity = Δn integer.")

print("\n5. THE FOCK JUMP OPERATOR — universal form, discrete kernel")
print(" 𝒥f(n) = Σ_z λ_z(n)[f(n+z) − f(n)] with z ∈ ℤ\{0}")
print(" bosons: z=±1 only (a, a†) -> nearest-rung jumps")
print(" pairs: z=±2 (a², a†²) -> squeezing, 2-photon processes")
print(" J-space: z=±1, n≤25 -> the workspace ladder, Δn=±1")
print(" -> EVERY quantum transition is 𝒥 with a discrete kernel. No exceptions.")
PY

[beanapologist]

No, not in main engine. as a separate example script in the example folder

Example: simulate the master equation in Fock space (diagonal part)

def jump_generator(N, gamma_up, gamma_down):
J = np.zeros((N, N))
for n in range(N):
if n > 0: J[n, n-1] += gamma_down * n # emission (z=-1)
if n < N-1: J[n, n+1] += gamma_up * (n+1) # absorption (z=+1)
np.fill_diagonal(J, -np.sum(J, axis=1))
return J