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

README.md

@@ -16,6 +16,22 @@ subject to continuous diffusion and discrete jump noise.
 - `f(Δ) = kΔ + gΔ²/(K²+Δ²)` — the nonlinear sink (linear + saturating)
 - `Δ* : Λ = f(Δ), f′(Δ*) > 0` — a stable equilibrium (basin centre)
 
+### The Jump Operator 𝒥
+
+Every Markov generator decomposes uniquely (**Lévy–Khintchine**) as `ℒ = b∂ₓ + ½σ²∂ₓ² + 𝒥` — drift, diffusion, jump, no fourth term. 𝒥 is the only term that evaluates *f at the destination* rather than derivatives at the origin.
+
+```
+𝒥f(x) = λ(x) ∫ [f(x+z) − f(x)] ν(dz|x)
+```
+
+`λ(x)` is the jump rate (`jump_rate`); `ν(dz|x)` is the jump kernel (`jump_size_dist`, default N(0,1)). Use `engine.jump_operator(f, x)` to evaluate 𝒥 numerically.
+
+**Confirmed properties** (verified by `TestJumpOperator`): annihilates constants · linear in `f` · quadratic identity `𝒥(x²) = λ` for N(0,1) kernel · zero under symmetric kernel for linear `f` · scales with rate · respects custom kernel.
+
+**Flat obstruction** — at a flat point `f(x) = e^{−1/x}` all derivatives vanish, so drift and diffusion go to zero. `𝒥f(0⁺) > 0` regardless. Local operators are blind to a flat wall; `𝒥` is not.
+
+**Crossing condition** — basin half-width `J_c = Δ_edge − Δ*`. Below `J_c`: `P(escape) ≈ 0`. Above: `P(escape) ≫ 0`. A cliff, not a slope.
+
 Use the `jump_diffusion_engine/` package to **analyse** stochastic systems and **steer** trajectories toward stable basins:
 
 | # | Action | Method | What it does |

examples/fock_ladder_generator.py

@@ -0,0 +1,167 @@
+"""
+fock_ladder_generator.py — Column generator for a Fock-space ladder
+
+Demonstrates five properties of the Markov generator matrix L for a
+truncated quantum harmonic oscillator / photon-number ladder:
+
+  dp/dt = L p,   with the column-generator convention  1^T L = 0
+
+The single-step transitions are standard bosonic ladder rates:
+  W_{n-1 <- n} = gd * n          (emission   / stimulated decay)
+  W_{n+1 <- n} = gu * (n+1)      (absorption / stimulated creation)
+
+Optional pair processes add dn=±2 jumps:
+  W_{n-2 <- n} = pair_dn * n*(n-1)       (two-photon loss)
+  W_{n+2 <- n} = pair_up * (n+1)*(n+2)   (two-photon gain)
+
+Demonstrations
+--------------
+1. COLUMN GENERATOR, dn=±1 only
+   - Columns sum to zero (1^T L = 0)
+   - Matrix is tridiagonal (bandwidth ≤ 1)
+   - Stationary distribution is geometric: p_n ∝ r^n, r = gu/gd
+
+2. BOUNDARY TERMS — truncation matters
+   - Lower wall n=0: emission rate gd*0 = 0  → natural reflecting wall
+   - Upper wall n=N-1: absorption would leave space, so it is dropped → reflecting
+   - n_bar converges to r/(1-r) as N → ∞; finite-N bias is truncation, not physics
+
+3. PAIR PROCESSES: dn=±2 widens the kernel support
+   - Bandwidth becomes 2; 1^T L = 0 still holds
+
+4. PARITY: pure dn=±2 splits Fock space into even/odd sectors
+   - Off-diagonal even↔odd block is exactly zero: parity is conserved
+   - Squeezing / two-photon processes cannot change photon-number parity
+
+5. DETAILED BALANCE with both dn=±1 and dn=±2 present
+   - Net current on each bond J_{a→b} = L[b,a]*p[a] − L[a,b]*p[b]
+   - Non-zero currents signal irreversibility when competing processes coexist
+
+Run:
+    python3 fock_ladder_generator.py
+"""
+
+import numpy as np
+
+
+def L_gen(N, gd, gu, pair_up=0.0, pair_dn=0.0):
+    """
+    Build the Markov generator matrix for a truncated Fock-space ladder.
+
+    Convention: L[n, m] = W_{n←m}  (off-diagonal, ≥ 0),
+                L[m, m] = −Σ_{n≠m} W_{n←m}  (diagonal, ≤ 0).
+
+    This is the *column-generator* form: dp/dt = L p, 1^T L = 0.
+
+    Parameters
+    ----------
+    N        : int   — number of Fock states  {0, 1, …, N-1}
+    gd       : float — single-photon loss rate (emission)
+    gu       : float — single-photon gain rate (absorption)
+    pair_up  : float — two-photon gain rate (dn = +2)
+    pair_dn  : float — two-photon loss rate (dn = −2)
+
+    Returns
+    -------
+    L : (N, N) ndarray
+    """
+    L = np.zeros((N, N))
+
+    def add(dest, src, rate):
+        """Add a transition src→dest with given rate (ignored if out of bounds)."""
+        if 0 <= dest < N and rate > 0:
+            L[dest, src] += rate
+            L[src, src]  -= rate
+
+    for n in range(N):
+        add(n - 1, n, gd * n)                  # emission:    dn = −1
+        add(n + 1, n, gu * (n + 1))             # absorption:  dn = +1
+        add(n - 2, n, pair_dn * n * (n - 1))   # pair loss:   dn = −2
+        add(n + 2, n, pair_up * (n + 1) * (n + 2))  # pair gain: dn = +2
+
+    return L
+
+
+def stationary(L):
+    """Return the normalised stationary distribution (eigenvector for eigenvalue 0)."""
+    w, v = np.linalg.eig(L)
+    i = np.argmin(np.abs(w))
+    p = np.real(v[:, i])
+    p /= p.sum()
+    return p
+
+
+def net_current(p, L, a, b):
+    """Net probability current on bond a↔b: J = L[b,a]*p[a] − L[a,b]*p[b]."""
+    return L[b, a] * p[a] - L[a, b] * p[b]
+
+
+# ---------------------------------------------------------------------------
+# 1. COLUMN GENERATOR, dn=±1 only
+# ---------------------------------------------------------------------------
+N = 8
+gd, gu = 1.0, 0.3
+
+L = L_gen(N, gd, gu)
+print("1. COLUMN GENERATOR, dn=+-1 only")
+print("   1^T L = ", np.round(L.sum(0), 12))
+
+idx = np.abs(np.subtract.outer(np.arange(N), np.arange(N)))
+bw_ok = np.all(idx[np.abs(L) > 1e-12] <= 1)
+print(f"   bandwidth: nonzero only on |n-m|<=1 -> {bw_ok}")
+
+p = stationary(L)
+r = gu / gd
+th = r ** np.arange(N)
+th /= th.sum()
+print(f"   stationary matches geometric r^n: {np.allclose(p, th)}")
+
+# ---------------------------------------------------------------------------
+# 2. BOUNDARY TERMS — truncation matters
+# ---------------------------------------------------------------------------
+print("\n2. BOUNDARY TERMS — the truncation matters")
+print("   at n=0: emission rate gd*0 = 0 (no leak below vacuum)  -> reflecting  OK")
+print(f"   at n=N-1: absorption gu*N = {gu*N:.2f} would leave the space; we DROP it.")
+print("   dropping = reflecting wall. This is why n_bar was slightly below r/(1-r):")
+
+nbar = np.sum(np.arange(N) * p)
+print(f"   n_bar(N={N}) = {nbar:.6f}   vs   r/(1-r) = {r/(1-r):.6f}")
+
+for NN in (8, 16, 32, 64):
+    LL = L_gen(NN, gd, gu)
+    pp = stationary(LL)
+    print(f"     N={NN:>3}: n_bar={np.sum(np.arange(NN)*pp):.9f}")
+
+print(f"   -> converges to {r/(1-r):.9f} as the wall recedes. Truncation, not physics.")
+
+# ---------------------------------------------------------------------------
+# 3. PAIR PROCESSES: dn=±2 widens the support of W
+# ---------------------------------------------------------------------------
+print("\n3. PAIR PROCESSES: dn=+-2 widens the support of W")
+L2 = L_gen(N, gd, gu, pair_up=0.02, pair_dn=0.05)
+bw = np.max(idx[np.abs(L2) > 1e-12])
+print(f"   bandwidth now = {bw} (was 1). Kernel support = {{+-1, +-2}}")
+print("   1^T L = ", np.round(L2.sum(0), 12))
+
+# ---------------------------------------------------------------------------
+# 4. PARITY: pure dn=±2 splits Fock space into even/odd sectors
+# ---------------------------------------------------------------------------
+print("\n4. PARITY: pure dn=+-2 splits Fock space into even/odd sectors")
+Lp = L_gen(N, 0, 0, pair_up=0.02, pair_dn=0.05)
+ev = np.arange(0, N, 2)
+od = np.arange(1, N, 2)
+cross = np.abs(Lp[np.ix_(od, ev)]).max()
+print(f"   coupling even->odd block: {cross:.3e}")
+print("   -> exactly zero: parity is CONSERVED. Two disconnected ladders.")
+print("      Squeezing/two-photon processes cannot change photon-number parity.")
+
+# ---------------------------------------------------------------------------
+# 5. DETAILED BALANCE with dn=±1 and dn=±2 present
+# ---------------------------------------------------------------------------
+print("\n5. DETAILED BALANCE with dn=+-2 present")
+p2 = stationary(L2)
+print("   net current on each bond (0 = reversible):")
+for n in range(4):
+    j1 = net_current(p2, L2, n, n + 1)
+    j2 = net_current(p2, L2, n, n + 2)
+    print(f"     {n}->{n+1}: {j1:+.3e}    {n}->{n+2}: {j2:+.3e}")

examples/master_equation.py

@@ -0,0 +1,158 @@
+"""
+master_equation.py — Pauli master equation via markov_generator
+
+Demonstrates that the generator matrix L built by JumpDiffusionEngine.markov_generator
+is the exact matrix encoding of the Pauli (Kolmogorov forward) master equation:
+
+    dp_n/dt = Σ_{m≠n} [ W_{n←m} p_m  −  W_{m←n} p_n ]
+                          ^^^^^^^^^^^^^    ^^^^^^^^^^^^^
+                              in              out
+
+In matrix form this is simply:
+
+    dp/dt = L · p
+
+where L[n, m] = W_{n←m}  (off-diagonal, ≥ 0)  and
+      L[n, n] = −Σ_{m≠n} W_{m←n} = −(rate_up[n] + rate_down[n])
+
+The diagonal is set by the *safer column-generator approach*:
+
+    L[n, n] = −(rate_up[n] + rate_down[n])
+
+so each column of L sums to exactly zero, which is the matrix-level statement
+of probability conservation: d/dt Σ_n p_n = 0.
+
+This script:
+  1. Builds L for the default engine (constant Λ, nonlinear sink f).
+  2. Starts from a narrow Gaussian initial distribution p(0).
+  3. Integrates dp/dt = L·p forward in time using scipy.integrate.solve_ivp.
+  4. Plots the evolving distribution and verifies that probability is conserved
+     at every output time.
+
+Run:
+    python3 master_equation.py
+"""
+
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from scipy.integrate import solve_ivp
+
+from jump_diffusion_engine import JumpDiffusionEngine
+
+
+# ---------------------------------------------------------------------------
+# 1. Set up the engine
+# ---------------------------------------------------------------------------
+LAMBDA_VAL = 0.5        # constant forcing Λ
+SIGMA      = 0.4        # diffusion strength
+K          = 0.8        # linear restoring rate
+G          = 0.5        # nonlinear gain
+K_SAT      = 2.0        # saturation scale
+
+eng = JumpDiffusionEngine(
+    lambda_func=lambda t: LAMBDA_VAL,
+    sigma=SIGMA,
+    jump_rate=0.0,       # pure diffusion + drift for clarity
+    dt=0.01,
+    seed=42,
+    k=K, g=G, K=K_SAT,
+)
+
+# ---------------------------------------------------------------------------
+# 2. Build the Markov generator matrix
+# ---------------------------------------------------------------------------
+X_RANGE = (-8.0, 12.0)
+N       = 300            # grid nodes  (matrix is N × N)
+
+out   = eng.markov_generator(lambda_val=LAMBDA_VAL, x_range=X_RANGE, n_points=N)
+L     = out['L']         # (N, N) generator — columns sum to zero
+x     = out['x']         # grid positions
+dx    = out['dx']
+
+# Verify the column-sum property (should be ~machine epsilon)
+max_col_err = np.abs(L.sum(axis=0)).max()
+print(f"Max |column sum| of L:  {max_col_err:.2e}  (should be < 1e-12)")
+
+# ---------------------------------------------------------------------------
+# 3. Initial distribution — narrow Gaussian centred at x = 2.0
+# ---------------------------------------------------------------------------
+x0_dist = 2.0            # initial mean
+sig0    = 0.4            # initial width
+
+p0 = np.exp(-0.5 * ((x - x0_dist) / sig0) ** 2)
+_trapz = np.trapezoid if hasattr(np, 'trapezoid') else np.trapz
+p0 /= _trapz(p0, x)   # normalise to a proper density
+
+# ---------------------------------------------------------------------------
+# 4. Integrate dp/dt = L · p forward in time
+# ---------------------------------------------------------------------------
+T_MAX    = 4.0
+T_REPORT = [0.0, 0.5, 1.0, 2.0, 4.0]   # snapshot times
+
+def rhs(t, p):
+    """dp/dt = L p"""
+    return L @ p
+
+sol = solve_ivp(
+    rhs,
+    t_span=(0.0, T_MAX),
+    y0=p0,
+    method='RK45',
+    t_eval=T_REPORT,
+    rtol=1e-6,
+    atol=1e-9,
+)
+
+assert sol.success, f"ODE solver failed: {sol.message}"
+
+# ---------------------------------------------------------------------------
+# 5. Verify probability conservation at every snapshot
+# ---------------------------------------------------------------------------
+print("\nProbability conservation check:")
+print(f"{'Time':>8s}   {'∫p dx':>12s}   {'max p':>10s}")
+for i, t in enumerate(sol.t):
+    p_t = sol.y[:, i]
+    norm = _trapz(p_t, x)
+    print(f"{t:8.2f}   {norm:12.6f}   {p_t.max():10.6f}")
+
+# ---------------------------------------------------------------------------
+# 6. Plot the evolving distribution
+# ---------------------------------------------------------------------------
+fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+
+# Left: probability distributions at each snapshot
+ax = axes[0]
+colours = plt.cm.viridis(np.linspace(0.1, 0.9, len(sol.t)))
+for i, t in enumerate(sol.t):
+    ax.plot(x, sol.y[:, i], color=colours[i], label=f"t = {t:.1f}")
+
+# Mark the stable fixed point
+fps = eng.find_fixed_points(LAMBDA_VAL)
+stable = [fp for fp in fps if fp['stable']]
+if stable:
+    xs = stable[0]['x_star']
+    ax.axvline(xs, color='crimson', linestyle='--', linewidth=1.2,
+               label=f"Δ* = {xs:.2f}")
+
+ax.set_xlabel("State Δ")
+ax.set_ylabel("Probability density p(Δ, t)")
+ax.set_title("Master equation: dp/dt = L · p")
+ax.legend(fontsize=8)
+ax.set_xlim(X_RANGE)
+ax.set_ylim(bottom=0)
+
+# Right: probability norm over time (conservation check)
+norms = [_trapz(sol.y[:, i], x) for i in range(len(sol.t))]
+ax2 = axes[1]
+ax2.plot(sol.t, norms, 'o-', color='steelblue')
+ax2.axhline(1.0, color='gray', linestyle='--', linewidth=0.8)
+ax2.set_xlabel("Time t")
+ax2.set_ylabel("∫ p(Δ, t) dΔ")
+ax2.set_title("Probability conservation (should stay = 1)")
+ax2.set_ylim(0.98, 1.02)
+
+plt.tight_layout()
+fig.savefig("/tmp/master_equation.png", dpi=120)
+print("\nFigure saved to /tmp/master_equation.png")