wai/wai-quantum-serve/Cargo.lock

@@ -589,7 +589,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.34"
+version = "0.3.35"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum/Cargo.lock

@@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.34"
+version = "0.3.35"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wai-quantum"
-version = "0.3.34"
+version = "0.3.35"
 edition = "2024"
 rust-version = "1.97"
 license = "Apache-2.0"
@@ -82,6 +82,9 @@ quantum_noisy = ["quantum_sv"]
 # Ground states of spin chains by two-site DMRG, held to exact diagonalisation
 # and to the free-fermion solution of the transverse-field Ising chain.
 quantum_dmrg = []
+# Time evolution of spin chains by the time-dependent variational principle,
+# held to exact evolution and to the free-fermion solution of a quench.
+quantum_tdvp = ["quantum_dmrg"]
 # Mixed states and the observables that distinguish quantum theory from a classical
 # account: density matrices, the partial trace, exact Pauli expectations, CHSH.
 quantum_info = ["quantum"]
@@ -121,7 +124,7 @@ quantum_phasor_meter = ["_signing"]
 
 # Everything.
 full = [
-    "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qio", "quantum_uf", "quantum_noisy", "quantum_dmrg", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel",
+    "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qio", "quantum_uf", "quantum_noisy", "quantum_dmrg", "quantum_tdvp", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel",
     "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec",
     "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom",
     "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata",

wai/wai-quantum/README.md

@@ -67,6 +67,7 @@ never creates information.
 | `quantum_qio` | **quantum-inspired optimisation**: Ising and Max-Cut problems (G-set reader, spin-glass generator), simulated bifurcation in its ballistic, discrete and heated forms (batched eight trials per pass, bit-identical to one at a time), simulated annealing, single-flip descent, an exact Gray-code referee up to 34 spins, and step-to-solution. Reaches the best-known cut on seven G-set instances |
 | `quantum_uf` | **union-find decoding in near-linear time** over the same matching graph: event-driven cluster growth on integer weights (a radix heap orders the fill events), merges, and peeling. Measured against exact matching on the same shots: within 0.92–1.18× the failures, up to 77× faster at `d = 11` |
 | `quantum_dmrg` | **ground states of spin chains by DMRG**: matrix-product operators for nearest-neighbour chains (XXZ Heisenberg in a field, transverse-field Ising), two-site sweeps with Lanczos and density-matrix truncation, entanglement entropies, a symmetric eigensolver, and two exact referees: full diagonalisation up to 22 sites, and the free-fermion solution of the Ising chain at any length. The critical Ising chain is reproduced to about `10⁻¹¹` per site at 200 sites |
+| `quantum_tdvp` | **time evolution of spin chains by TDVP**: two-site sweeps forward and back, Lanczos exponentials that halve a step they cannot converge, bonds split by singular value decomposition, and every basis completed with null vectors so a product state's small bonds do not hold the evolution back. Measures site and bond observables, entanglement at every bond, energy and the Loschmidt echo. Held to two exact referees: Krylov evolution of all amplitudes up to 20 sites, and the free-fermion solution of an Ising quench at any length. On a 100-site critical quench at bond 64 it stays within `10⁻⁹` of the exact answer until the entanglement fills the bond, and past that its error follows the discarded weight it reports |
 | `quantum_noisy` | **noisy circuits**: quantum trajectories on the state vector (depolarizing, dephasing, amplitude damping, readout), an exact density-matrix referee up to 12 qubits that reuses the same kernels, linear cross-entropy fidelity (sampled, exact and the circuit's own), and the digital error model. Measures where that model is off: below the true fidelity by a margin that grows with the error rate, not with depth |
 | `quantum_ldpc` | **qLDPC memory at circuit level**: bivariate-bicycle codes from their two polynomials (`k` by rank, logical bases), the depth-8 syndrome cycle with circuit-level noise, its detector error model, and BP-OSD decoding (min-sum with adaptive scaling, ordered statistics with a combination sweep). The error model has the same columns as the source's own decoding matrices, and on identical syndromes the decoder makes the source's decisions shot for shot |
 | `quantum_arch` | **architecture-level estimation**: yoked surface-code storage (row and grid parity checks over patches, hot and cold), magic-state cultivation read off the published simulation curves, CCZ factories, and a clock set by whichever is slowest of lattice surgery, the magic-state supply and the control loop. Assembles whole-machine estimates; its storage search matches the source's footprint tool layout for layout |
@@ -197,6 +198,28 @@ At 200 sites the critical Ising chain comes out within about `10⁻¹¹` per sit
 of its exact energy. Runs repeat bit for bit at any thread count, natively and
 as WebAssembly.
 
+## Time evolution
+
+```sh
+wai-quantum tdvp -n 100 --bond 64 --dt 0.1 --t 4             # critical Ising quench, against free fermions
+wai-quantum tdvp --chain heisenberg -n 100 --bond 64 --t 3    # the Néel state melting, entanglement included
+```
+
+A quantum simulation of dynamics is checkable when a classical method can
+follow it. `quantum_tdvp` evolves a matrix-product state by the
+time-dependent variational principle and reports what it gave up: the
+discarded weight, the energy drift and the entanglement at every bond. It is
+held to exact evolution of every amplitude on short chains, and to the
+free-fermion solution of the Ising quench at any length, including the
+Loschmidt echo.
+
+On the 100-site critical Ising quench at bond 64, every site's `σˣ` stays
+within `10⁻⁹` of the exact answer through `t = 3`, with the Loschmidt rate
+agreeing to seven digits. When the entanglement outgrows the bond, the error
+rises with the discarded weight, and the run says so. A chain small enough to
+fit the bond evolves exactly, to rounding, at any time step. Runs repeat bit
+for bit at any thread count, natively and as WebAssembly.
+
 ## Noisy circuits
 
 ```sh

wai/wai-quantum/examples/tdvp_quench.rs

@@ -0,0 +1,54 @@
+//! A quench of the transverse-field Ising chain by TDVP, held at every
+//! report to the exact free-fermion solution: the chain starts with every
+//! σˣ = +1 and evolves under `−J Σ σᶻσᶻ − h Σ σˣ`.
+//!
+//! ```text
+//! cargo run --release --example tdvp_quench --features quantum_tdvp -- N BOND DT T [h] [cutoff] [threads]
+//! ```
+use wai_quantum::quantum_dmrg::{Chain, Mps, Op, Scalar, C};
+use wai_quantum::quantum_tdvp::{ising_quench, Evolution, TdvpConfig};
+
+fn main() {
+    let a: Vec<String> = std::env::args().collect();
+    let n: usize = a[1].parse().unwrap();
+    let bond: usize = a[2].parse().unwrap();
+    let dt: f64 = a[3].parse().unwrap();
+    let t_end: f64 = a[4].parse().unwrap();
+    let h: f64 = a.get(5).map_or(1.0, |x| x.parse().unwrap());
+    let cutoff: f64 = a.get(6).map_or(1e-20, |x| x.parse().unwrap());
+    let threads: usize = a.get(7).map_or(TdvpConfig::default().threads, |x| x.parse().unwrap());
+    let chain = Chain::ising(n, 1.0, h);
+    let plus = Mps::product(&vec![[C { re: 1.0, im: 0.0 }, C { re: 1.0, im: 0.0 }]; n]);
+    let cfg = TdvpConfig { dt, max_bond: bond, cutoff, threads, ..TdvpConfig::default() };
+    let mut ev = Evolution::new(&chain, &plus, cfg);
+    let e0 = ev.energy();
+    let steps = (t_end / dt).round() as usize;
+    let every = ((0.5 / dt).round() as usize).max(1);
+    let t0 = std::time::Instant::now();
+    let mut discarded = 0.0;
+    println!("ising n={n} h={h} bond={bond} dt={dt} cutoff={cutoff:e}");
+    println!("    t   max|Δσˣ|   max|Δσᶻσᶻ|  σˣ(mid)        exact          rate      exact rate  S(mid)  bond  |ΔE|     discarded  seconds");
+    for s in 1..=steps {
+        discarded += ev.step().discarded;
+        if s % every == 0 || s == steps {
+            let t = ev.time();
+            let ff = ising_quench(n, 1.0, h, t);
+            let x = ev.expect(Op::x());
+            let zz = ev.bond_expect(Op::z(), Op::z());
+            let dx = x.iter().zip(&ff.x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max);
+            let dzz = zz.iter().zip(&ff.zz).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max);
+            let echo = ev.overlap(&plus).norm2();
+            let rate = -echo.ln() / n as f64;
+            let exact_rate = -ff.echo.ln() / n as f64;
+            let s_mid = ev.entropies()[n / 2 - 1];
+            println!(
+                "{t:5.2}  {dx:.2e}   {dzz:.2e}    {:+.10}  {:+.10}  {rate:.8}  {exact_rate:.8}  {s_mid:.4}  {:4}  {:.1e}  {discarded:.1e}  {:.1}",
+                x[n / 2],
+                ff.x[n / 2],
+                ev.max_bond(),
+                (ev.energy() - e0).abs(),
+                t0.elapsed().as_secs_f64()
+            );
+        }
+    }
+}

wai/wai-quantum/src/lib.rs

@@ -17,7 +17,7 @@
 #[cfg(feature = "_signing")]
 mod merkle;
 
-#[cfg(any(feature = "quantum_mps", feature = "quantum_bptn", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_sv"))]
+#[cfg(any(feature = "quantum_mps", feature = "quantum_bptn", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_sv", feature = "quantum_dmrg"))]
 mod linalg;
 
 // Gates as matrices, shared by the state vector and the tensor networks.
@@ -38,7 +38,8 @@ pub mod gates;
     feature = "quantum_sv",
     feature = "quantum_arch",
     feature = "quantum_qio",
-    feature = "quantum_dmrg"
+    feature = "quantum_dmrg",
+    feature = "quantum_tdvp"
 ))]
 mod repro;
 
@@ -123,6 +124,9 @@ pub mod quantum_noisy;
 #[cfg(feature = "quantum_dmrg")]
 pub mod quantum_dmrg;
 
+#[cfg(feature = "quantum_tdvp")]
+pub mod quantum_tdvp;
+
 #[cfg(feature = "quantum_frame")]
 pub mod quantum_match;
 

wai/wai-quantum/src/linalg.rs

@@ -8,7 +8,7 @@
 // Complex f64 (deterministic: only + - * / sqrt, no FMA/transcendentals)
 // ===========================================================================
 
-#[derive(Clone, Copy, Debug)]
+#[derive(Clone, Copy, Debug, PartialEq)]
 pub struct C {
     pub re: f64,
     pub im: f64,
@@ -148,7 +148,7 @@ pub(crate) fn jacobi_svd(mut m: Vec<Vec<C>>, rows: usize, cols: usize) -> (Vec<V
 /// rotates nothing. The original keeps its looser stopping rule (columns
 /// orthogonal to about 1e-7 relative) because the MPS backend's results are
 /// pinned to it.
-#[cfg(feature = "quantum_bptn")]
+#[cfg(any(feature = "quantum_bptn", feature = "quantum_tdvp"))]
 pub(crate) fn jacobi_svd_strict(mut m: Vec<Vec<C>>, rows: usize, cols: usize) -> (Vec<Vec<C>>, Vec<f64>, Vec<Vec<C>>) {
     // V starts as identity (cols × cols), stored by columns.
     let mut v: Vec<Vec<C>> = (0..cols)

wai/wai-quantum/src/repro.rs

@@ -15,7 +15,8 @@
     feature = "quantum_phasor",
     feature = "quantum_kernel",
     feature = "quantum_tn",
-    feature = "quantum_sv"
+    feature = "quantum_sv",
+    feature = "quantum_tdvp"
 ))]
 mod trig;
 #[cfg(any(
@@ -26,7 +27,8 @@ mod trig;
     feature = "quantum_phasor",
     feature = "quantum_kernel",
     feature = "quantum_tn",
-    feature = "quantum_sv"
+    feature = "quantum_sv",
+    feature = "quantum_tdvp"
 ))]
 pub use trig::sin_cos;
 

wai/wai-quantum/tests/bit_identity.rs

@@ -25,7 +25,8 @@
     feature = "quantum_qio",
     feature = "quantum_uf",
     feature = "quantum_noisy",
-    feature = "quantum_dmrg"
+    feature = "quantum_dmrg",
+    feature = "quantum_tdvp"
 ))]
 
 use wai_quantum::*;
@@ -474,3 +475,36 @@ fn dmrg_ground_states() {
     v.extend(vecs.iter().map(|x| x.abs()));
     assert_eq!(digest(&v), 0x3b28_8bdd_be3f_48a2);
 }
+
+#[test]
+fn tdvp_time_evolution() {
+    use quantum_dmrg::{Chain, Mps, Op, Scalar, C};
+    use quantum_tdvp::*;
+    let mut v = Vec::new();
+    // A truncated Heisenberg quench from the Néel state, and the Ising quench
+    // with its free-fermion referee.
+    let n = 10;
+    let neel = Mps::basis(&(0..n).map(|i| i % 2 == 1).collect::<Vec<_>>());
+    let mut ev = Evolution::new(&Chain::heisenberg(n, 1.0, 0.7, 0.1), &neel, TdvpConfig { dt: 0.1, max_bond: 6, threads: 1, ..TdvpConfig::default() });
+    for _ in 0..6 {
+        v.push(ev.step().discarded);
+    }
+    v.extend(ev.expect(Op::z()));
+    v.extend(ev.bond_expect(Op::sz(), Op::sz()));
+    v.extend(ev.entropies());
+    v.push(ev.energy());
+    let o = ev.overlap(&neel);
+    v.extend([o.re, o.im]);
+    let plus = Mps::product(&vec![[C::ONE, C::ONE]; 8]);
+    let mut ev = Evolution::new(&Chain::ising(8, 1.0, 0.9), &plus, TdvpConfig { dt: 0.2, max_bond: 16, cutoff: 0.0, threads: 1, ..TdvpConfig::default() });
+    for _ in 0..3 {
+        ev.step();
+    }
+    v.extend(ev.expect(Op::x()));
+    let ff = ising_quench(8, 1.0, 0.9, ev.time());
+    v.extend(ff.x);
+    v.extend(ff.zz);
+    v.push(ff.echo);
+    assert_eq!(digest(&v), 0xe1b4_fd28_af23_50cd);
+}
+

wai/wai-rs/Cargo.lock

@@ -3900,7 +3900,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.34"
+version = "0.3.35"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-rs/Cargo.toml

@@ -368,7 +368,7 @@ quantum_arch = ["dep:wai-quantum", "wai-quantum/quantum_arch"]
 # qLDPC memory at circuit level: bivariate-bicycle codes, BP-OSD.
 quantum_ldpc = ["quantum_frame", "wai-quantum/quantum_ldpc"]
 # Quantum-inspired optimisation: simulated bifurcation, annealing, exact referee.
-quantum_qio = ["dep:wai-quantum", "wai-quantum/quantum_qio", "wai-quantum/quantum_dmrg"]
+quantum_qio = ["dep:wai-quantum", "wai-quantum/quantum_qio", "wai-quantum/quantum_dmrg", "wai-quantum/quantum_tdvp"]
 quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"]
 # Variational quantum machine learning (extensions/quantum-sim § VML,
 # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot

wai/wai-rs/src/bin/wai_quantum.rs

@@ -167,7 +167,9 @@ fn usage() -> ! {
          wai-quantum gbs … --check SAMPLES.txt        (observed patterns vs exact marginals)\n\n\
          SPIN CHAINS (DMRG ground states, exact referees):\n  \
          wai-quantum dmrg [--chain ising|heisenberg] [-n N] [--j J] [--jz JZ] [--h H] [--bond D] [--sweeps S]\n    \
-         [--cutoff C] [--exact] [--json]\n\n\
+         [--cutoff C] [--exact] [--json]\n  \
+         wai-quantum tdvp [--chain ising|heisenberg] [-n N] [--j J] [--jz JZ] [--h H] [--start plus|neel|up]\n    \
+         [--bond D] [--dt DT] [--t T] [--every K] [--cutoff C] [--json]   (time evolution by TDVP)\n\n\
          OPTIMISATION (simulated bifurcation, annealing, exact):\n  \
          wai-quantum qio (--gset FILE | --sk N [--instance S]) [--method dsb|bsb|hdsb|hbsb|sa|exact]\n    \
          [--steps N] [--trials T] [--dt D] [--c1 C] [--gamma G] [--beta0 B --beta1 B] [--descend]\n    \
@@ -1030,6 +1032,90 @@ fn main() {
             }
         }
 
+        // Time evolution of a spin chain by TDVP from a product state. The
+        // Ising chain started with every σˣ = +1 is held at each report to
+        // its free-fermion solution.
+        "tdvp" => {
+            use wai::quantum_dmrg::{Chain, Mps, Op, Scalar, C};
+            use wai::quantum_tdvp::*;
+            let n: usize = num(args, "--sites", "-n", 40);
+            if n < 2 {
+                die("-n: at least 2 sites");
+            }
+            let kind = opt(args, "--chain", "").unwrap_or_else(|| "ising".into());
+            let j: f64 = num(args, "--j", "", 1.0);
+            let h: f64 = num(args, "--h", "", if kind == "heisenberg" { 0.0 } else { 1.0 });
+            let chain = match kind.as_str() {
+                "ising" => Chain::ising(n, j, h),
+                "heisenberg" => Chain::heisenberg(n, j, num(args, "--jz", "", j), h),
+                _ => die("--chain: ising or heisenberg"),
+            };
+            let start_kind = opt(args, "--start", "").unwrap_or_else(|| if kind == "ising" { "plus".into() } else { "neel".into() });
+            let start = match start_kind.as_str() {
+                "plus" => Mps::product(&vec![[C::ONE, C::ONE]; n]),
+                "neel" => Mps::basis(&(0..n).map(|i| i % 2 == 1).collect::<Vec<_>>()),
+                "up" => Mps::basis(&vec![false; n]),
+                _ => die("--start: plus, neel or up"),
+            };
+            let dt: f64 = num(args, "--dt", "", 0.05);
+            let t_end: f64 = num(args, "--t", "", 2.0);
+            if dt.is_nan() || dt <= 0.0 || t_end.is_nan() || t_end < 0.0 {
+                die("--dt must be positive and --t not negative");
+            }
+            let cfg = TdvpConfig { dt, max_bond: num(args, "--bond", "", 32), cutoff: num(args, "--cutoff", "", 1e-20), ..TdvpConfig::default() };
+            let steps = (t_end / dt).round() as usize;
+            let every: usize = num(args, "--every", "", ((0.5 / dt).round() as usize).max(1));
+            // The observable at the middle site: σˣ for the Ising chain, σᶻ
+            // for the Heisenberg chain.
+            let (obs, obs_name) = if kind == "ising" { (Op::x(), "x") } else { (Op::z(), "z") };
+            let referee = kind == "ising" && start_kind == "plus";
+            let mut ev = Evolution::new(&chain, &start, cfg);
+            let e0 = ev.energy();
+            let mut discarded = 0.0;
+            let mut rows: Vec<String> = Vec::new();
+            let json = flag(args, "--json");
+            if !json {
+                println!("{kind} chain, {n} sites, start {start_kind}, bond ≤ {}, dt {dt}", cfg.max_bond);
+                print!("     t   {obs_name}(mid)         ");
+                if referee {
+                    print!("exact          max error  rate        exact rate  ");
+                }
+                println!("S(mid)   bond  |ΔE|     discarded");
+            }
+            for s in 1..=steps {
+                discarded += ev.step().discarded;
+                if s % every != 0 && s != steps {
+                    continue;
+                }
+                let t = ev.time();
+                let o = ev.expect(obs);
+                let s_mid = ev.entropies()[n / 2 - 1];
+                let de = (ev.energy() - e0).abs();
+                let rate = -ev.overlap(&start).norm2().ln() / n as f64;
+                let ff = referee.then(|| ising_quench(n, j, h, t));
+                let max_err = ff.as_ref().map(|f| o.iter().zip(&f.x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max));
+                if json {
+                    rows.push(format!(
+                        "{{\"bond\":{},\"discarded\":{discarded},\"energy_drift\":{de},\"entropy\":{s_mid},\"exact\":{},\"exact_rate\":{},\"max_error\":{},\"mid\":{},\"rate\":{rate},\"t\":{t}}}",
+                        ev.max_bond(),
+                        ff.as_ref().map_or("null".into(), |f| f.x[n / 2].to_string()),
+                        ff.as_ref().map_or("null".into(), |f| (-f.echo.ln() / n as f64).to_string()),
+                        max_err.map_or("null".into(), |e| e.to_string()),
+                        o[n / 2]
+                    ));
+                } else {
+                    print!("{t:6.2}  {:+.10}  ", o[n / 2]);
+                    if let (Some(f), Some(e)) = (&ff, max_err) {
+                        print!("{:+.10}  {e:.2e}   {rate:.8}  {:.8}  ", f.x[n / 2], -f.echo.ln() / n as f64);
+                    }
+                    println!("{s_mid:.4}  {:4}  {de:.1e}  {discarded:.1e}", ev.max_bond());
+                }
+            }
+            if json {
+                println!("{{\"chain\":\"{kind}\",\"dt\":{dt},\"rows\":[{}],\"sites\":{n},\"start\":\"{start_kind}\"}}", rows.join(","));
+            }
+        }
+
         // Quantum-inspired optimisation: Max-Cut (G-set text) or a spin
         // glass, by simulated bifurcation, annealing, or exact enumeration.
         // Max-Cut reports cuts; a spin glass reports energies.

wai/wai-rs/src/lib.rs

@@ -663,6 +663,9 @@ pub use wai_quantum::quantum_qio;
 #[cfg(feature = "quantum_qio")]
 pub use wai_quantum::quantum_dmrg;
 
+#[cfg(feature = "quantum_qio")]
+pub use wai_quantum::quantum_tdvp;
+
 #[cfg(any(feature = "quantum_tn", feature = "quantum_sv"))]
 pub use wai_quantum::gates;
 #[cfg(feature = "quantum_qir")]

wai/wai-web/Cargo.lock

@@ -734,7 +734,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.34"
+version = "0.3.35"
 dependencies = [
  "blake3",
  "ed25519-dalek",