Gibbs State Preparation¶

Category C2 · Module oracq.algorithms.input_model.density · Stage V1

Overview¶

Given a block encoding of a Hamiltonian \(H\) and an inverse temperature \(\beta \ge 0\), prepares an approximate purification of the Gibbs state \(\rho = e^{-\beta H}/Z\) (\(Z = \operatorname{Tr} e^{-\beta H}\)). The implementation takes the QSVT purification route (Chowdhury–Somma 2017, van Apeldoorn–Gilyén 2019, Gilyén et al. 2019, arXiv:1806.01838): first prepare the purification of the maximally mixed state on system and environment (\(n\) Bell pairs, see Purification access), then apply to system a QSVT block encoding whose zero-signal block is proportional to \(g(H/\alpha)\), where on the spectral variable \(x = \lambda/\alpha \in [-1,1]\)

\[ g(x) = \exp\!\bigl(-\tfrac{\beta\alpha}{2}(x+1)\bigr) \in (0, 1], \qquad c = \beta\alpha/2 . \]

\(g\) is split by parity into the two branches \(e^{-c}\cosh(cx)\) and \(-e^{-c}\sinh(cx)\), each approximated by a modified Bessel truncation; after post-selecting signal == 0, the reduced density matrix of system is proportional to \(g(H/\alpha)^2 = e^{-\beta H}\) (up to normalization).

Interface and input model¶

gibbs_purification(hamiltonian, beta, *, error=0.01)

API entry points: gibbs_purification

hamiltonian: a BlockEncoding; by convention the spectrum of \(H\) lies in \([-\alpha, \alpha]\) (input model BE + DM; \(\alpha\) = be_alpha).
beta: the inverse temperature \(\beta \ge 0\); at \(\beta = 0\) it degenerates to the maximally mixed purification (qsp_degree = 0, gibbs_scale = 1.0).
error: the uniform polynomial approximation error, which must lie in \((0, 1)\); each branch’s truncation tail is controlled to \(\le\) error/8.

Returns an ApproximatePurification (oracle_kind = "approximate_purification", operation signature ("system", "environment", "signal")). Module attributes:

Attribute	Meaning
`algorithm`	`"gibbs_purification"`
`beta` / `error`	echo of the call parameters
`qsp_degree`	the larger of the even/odd branch truncation degrees
`gibbs_scale`	the zero-signal block normalization \(2s\) (independent of the partition function; does not affect the reduced state)
`success_condition`	`"signal == 0"`
`width` / `environment_width` / `signal_qubits`	widths of the three registers

The classical reference gibbs_state(hamiltonian, beta) (small-matrix \(e^{-\beta H}/\operatorname{Tr}\)) and trace_distance serve the witnesses and downstream reuse.

Implementation notes¶

The generation chain is Bell-pair preparation → even/odd branch truncation → phase synthesis → LCU addition. Both branches are convex; the endpoint-matched imaginary-part completion always satisfies the unit-disk constraint (\(f^2(x)\) does not exceed the chord joining the endpoints), so phase synthesis is always feasible; each branch, after its phases are synthesized via imaginary-part completion, has its real part extracted with \((U_\Phi + U_{-\Phi})/2\), and the two branches are then added via linear_combination. The scaling is \(s = 1.5\max(\lVert g_{\mathrm{even}}\rVert_\infty, \lVert g_{\mathrm{odd}}\rVert_\infty, 10^{-3})\), with the two branches’ combined uniform error at most error/4.

Register layout: system(n) | environment(n) | signal(gibbs_be.signal_qubits), where environment carries the other half of the Bell pairs. Applicability boundary: when \(\beta\alpha\) is too large the polynomial degree exceeds the synthesis cap (40), and a ValidationError is raised at generation time (suggesting a smaller \(\beta\) or first shrinking the spectral scale of \(H\)), never silently degrading; the input must be a block encoding, not a sparse oracle.

Validation approach¶

Category C2 (approximately continuous semantics; the acceptance criteria are in ../development/validation-plan.md §2): the post-selected reduced density matrix must approximate the classical reference Gibbs state within tolerance. Three layers of evidence:

Structure: tests/core/test_density.py:GibbsTests.test_beta_zero_is_maximally_mixed (attributes of the \(\beta=0\) degenerate path and a cross-check against \(I/2\)) and test_invalid_inputs_fail_at_generation (\(\beta < 0\), error equal to 0 or \(\ge 1\), a hamiltonian that is not a BlockEncoding, etc.).
Numerical: GibbsTests.test_diagonal_hamiltonian_gibbs and test_non_diagonal_hamiltonian_gibbs — for 2×2 diagonal / non-diagonal Hamiltonians (\(\beta = 0.6 / 0.8\), error = 0.05), the signal == 0 branch is taken on the reference simulator and, after the partial trace and renormalization, cross-checked against the gibbs_state classical reference with trace distance < 0.1.
Binding: this algorithm has no separate binding witness (the input already requires a concrete BE).

Convergence witnesses (added in V1): test_error_convergence_decreases fixes \(\beta = 0.8\) and sweeps error ∈ {0.4, 0.2, 0.1}, asserting that the trace distance is monotonically non-increasing across levels and that each level is ≤ error — measured 5.6e-4 / 8.7e-5 / 8.7e-5 (error 0.2 and 0.1 land on the same truncation degree and give equal distances, hence the assertion is non-increasing rather than strictly decreasing); test_error_bound_uniform_in_beta fixes error = 0.1 and sweeps \(\beta\) ∈ {0.2, 0.5, 1.0}, with each point’s trace distance ≤ error — measured 2.7e-6 / 1.5e-5 / 1.9e-4, growing with \(c = \beta\alpha/2\) but always far below error.

Known gaps and planned stages¶

No known gaps; the stage V1 witnesses are complete (diagonal/non-diagonal cross-checks + a monotonically non-increasing sweep over three error levels + β-grid consistency).

Numerical validation¶

Paper-grade numerical experiments are in tests/verification/verify_hamiltonian.py (real-backend execution, no mock substitutes). Experiment design: the non-diagonal Hamiltonian \(H = \begin{pmatrix} 0.5 & 0.2 \\ 0.2 & -0.3 \end{pmatrix}\) (\(\alpha = 0.7\), \(\beta = 1.2\)) and the diagonal Hamiltonian \(H = \operatorname{diag}(0.8, -0.4)\) (\(\beta = 0.6\)), with error = 0.02; plus a \(\beta = 0\) degenerate case. The purification program is executed on the reference / rir-pysparq / originir-ext paths, the signal == 0 branch is taken, the partial trace is computed with an independent numpy implementation and renormalized, the trace distance against the classical Gibbs state \(e^{-\beta H}/Z\) computed by scipy.linalg.expm is compared (numpy.linalg.eigvalsh), and the post-selection success probability is reported.

Case	Scale	Paths	Metric	Value
`gibbs-trace-distance-nondiag`	13 qubits, \(\beta = 1.2\)	reference, rir-pysparq, originir-ext	trace distance	5.10e-5
ditto	—	—	success probability	0.0953
`gibbs-trace-distance-diag`	diagonal \(H\), \(\beta = 0.6\)	reference, rir-pysparq	trace distance	1.28e-5
ditto	—	—	success probability	0.0992
`gibbs-beta0-maximally-mixed`	\(\beta = 0\) degenerate path	reference, rir-pysparq, originir-ext	trace distance	0.0 (exactly \(I/2\))

The trace distances of all three cases are far below the 3×error = 0.06 criterion (highest measured 5.1e-5), consistent with the convergence witnesses recorded on this page (the trace distance is monotonically non-increasing in error and far below error).

Reproduction command:

PYTHONPATH=src /path/to/backend/python tests/verification/verify_hamiltonian.py

Artifacts: out/verification/hamiltonian.json.