Complex interpolation of matrices from multi-manifold learning can be used to enhance the analysis of ergodicity onset in disordered quantum systems simulated on digital quantum processors.

Given noisy time-series or shadow-tomography data collected from a digital quantum processor executing a disordered spin-chain circuit (e.g., a kicked Ising or Heisenberg model with tunable disorder strength W), constructing a continuous family of effective density matrices via complex (geodesic/log-Euclidean or Riemannian) interpolation across multiple learned data manifolds (one per disorder realization or noise stratum) will yield an ergodicity-onset indicator — e.g., interpolated level-spacing ratio ⟨r⟩, interpolated bipartite entanglement entropy, or an interpolated out-of-time-order correlator (OTOC) — whose estimated critical disorder strength W_c(interp) differs from the estimate obtained via direct (non-interpolated) hardware averaging by ≥30% smaller bias relative to a noiseless statevector-simulator ground truth, at fixed shot budget and fixed hardware noise level (measured via randomized benchmarking, at least 1×10⁻³ per-gate error).

Minimum viable test (MVT): N=10 qubits, single disorder-driven Floquet/kicked-Ising circuit family, 5 disorder strengths W spanning the expected crossover, 3 manifolds (three independent noise-mitigation strata: raw, ZNE-2x, ZNE-3x fold), one real backend (e.g., IBM Heron r2 or similar ≥99% median 2Q gate fidelity device) plus one noiseless statevector simulator as ground truth. Estimate ⟨r⟩ or half-chain entanglement entropy at each W via 4,096 shots per circuit, 20 disorder realizations per W. Compare interpolated vs. direct-averaged estimator bias against simulator truth using bootstrap resampling (1,000 resamples) for confidence intervals.

# Pseudocode
for W in disorder_strengths:
    gt[W] = exact_diagonalization(N, W, n_realizations=200)  # ground truth

manifolds = ['raw', 'zne_2x', 'alt_transpile']
hw_data = {}
for backend in [backend_A, backend_B]:
    calibrate_and_run_RB(backend)
    for W in disorder_strengths:
        for manifold in manifolds:
            circuits = build_disordered_circuits(N, W, manifold, n_realizations=20)
            counts = execute_on_hardware(backend, circuits, shots=4096)
            rho_est = reconstruct_density_matrices(counts)  # via shadow tomography
            hw_data[(backend, W, manifold)] = rho_est

    for W in disorder_strengths:
        embeddings = {m: manifold_embed(hw_data[(backend, W, m)]) for m in manifolds}
        rho_interp = riemannian_interpolate(embeddings, metric='log_euclidean')
        est_interp[W] = compute_ergodicity_indicator(rho_interp)   # <r>, entropy, OTOC
        est_baseline[W] = compute_ergodicity_indicator(
            weighted_average([hw_data[(backend,W,m)] for m in manifolds]))

    Wc_interp = fit_sigmoid_crossover(disorder_strengths, est_interp)
    Wc_baseline = fit_sigmoid_crossover(disorder_strengths, est_baseline)
    Wc_gt = fit_sigmoid_crossover(disorder_strengths, gt)

    bias_interp = abs(Wc_interp - Wc_gt)
    bias_baseline = abs(Wc_baseline - Wc_gt)
    bootstrap_compare(bias_interp, bias_baseline, n_resamples=1000)