Quantum Gibbs Sampler

• Quantum Gibbs samplers are quantum algorithms that prepare thermal Gibbs states of Hamiltonians using open-system dynamics and circuit-based methods.
• They leverage techniques like Lindbladian evolution, phase estimation, and amplitude amplification to control error and mixing time.
• Applications include estimating thermodynamic quantities, benchmarking quantum hardware, and demonstrating quantum computational advantage.

A quantum Gibbs sampler is a quantum algorithm or open-system dynamics for preparing the Gibbs state of a Hamiltonian $H$ at inverse temperature $\beta$, typically written as $\rho_\beta=e^{-\beta H}/Z$ with $Z=\mathrm{Tr}[e^{-\beta H}]$. For a many-body Hamiltonian $H=\sum_i h_i$ on Hilbert-space dimension $D$, the operational task is to output a state $\rho$ such that $\|\rho-\rho_\beta\|_1\le \epsilon$ while controlling implementation cost, locality, and mixing time (Zhang et al., 2023). In the dissipative formulation, a quantum Gibbs sampler is a local, primitive, reversible Lindblad generator whose unique steady state is the Gibbs state; in more recent noncommutative formulations, the sampler is engineered to satisfy exact Kubo–Martin–Schwinger detailed balance for arbitrary Hamiltonians (Kastoryano et al., 2014, Chen et al., 2023).

1. Formal setting and target objects

The central object is the thermal state

$$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$

which describes equilibrium at temperature $T=\beta^{-1}$. In algorithmic settings, one asks either for a mixed-state preparation procedure, for a sampler that outputs measurement results distributed according to $\beta$0, or for a purified preparation of the thermofield-double state

$\beta$1

whose partial trace returns $\beta$2 (1603.02940, Leng et al., 24 Apr 2026).

Two notions of “sampler” coexist. In the semigroup viewpoint, one studies $\beta$3 for a Lindbladian $\beta$4 with $\beta$5, and convergence is governed by spectral gap, log-Sobolev behavior, or related mixing parameters. In the circuit viewpoint, one prepares $\beta$6 or a purification directly using Hamiltonian simulation, phase estimation, linear combinations of unitaries, adiabatic continuation, or measurement-based updates. This suggests that “quantum Gibbs sampler” is not a single algorithmic template but a family of constructions sharing a common equilibrium target.

The distinction between state preparation and sample generation is not merely terminological. Some works target the density operator itself, some target computational-basis samples from a Gibbs state, and some exploit the sampler to estimate thermodynamic quantities such as $\beta$7 or the free energy. That breadth is already visible in early phase-estimation-based algorithms, dissipative Lindbladians, and recent continuous-variable constructions (1603.02940, Becker et al., 16 Apr 2026).

2. Detailed balance, Lindbladians, and equilibrium structure

For open-system samplers, the basic object is a Gorini–Kossakowski–Sudarshan–Lindblad generator. In the Heisenberg picture, a general KMS-detailed-balanced form can be written as

$\beta$8

with the KMS condition $\beta$9, equivalently

$\rho_\beta=e^{-\beta H}/Z$0

A structural characterization states that $\rho_\beta=e^{-\beta H}/Z$1 satisfies $\rho_\beta=e^{-\beta H}/Z$2-KMS detailed balance iff it admits such a representation with a coherent term $\rho_\beta=e^{-\beta H}/Z$3 that commutes with $\rho_\beta=e^{-\beta H}/Z$4 (Ding et al., 2024).

This framework generalizes the older commuting-Hamiltonian theory. For commuting local Hamiltonians, the Davies generator and the heat-bath generator are both local, primitive, reversible, frustration-free Gibbs samplers. Their convergence is governed by the spectral gap

$\rho_\beta=e^{-\beta H}/Z$5

and a central equivalence states that the gap is independent of system size iff the Gibbs state satisfies a strong clustering property. In one dimension at any temperature, and in any dimension at high enough temperature, the samplers are gapped, so the Gibbs state can be prepared in polynomial time on a quantum computer (Kastoryano et al., 2014).

For arbitrary noncommutative Hamiltonians, Chen, Kastoryano, and Gilyén construct an exact detailed-balanced Lindbladian

$\rho_\beta=e^{-\beta H}/Z$6

where $\rho_\beta=e^{-\beta H}/Z$7 are frequency-filtered Heisenberg evolutions of local scrambling operators and $\rho_\beta=e^{-\beta H}/Z$8 is a coherent correction chosen so that exact quantum detailed balance holds. For lattice Hamiltonians, Lieb–Robinson bounds imply that the dressed jump operators are quasi-local with radius $\rho_\beta=e^{-\beta H}/Z$9, so the algorithm inherits a local-update interpretation even for noncommuting systems (Chen et al., 2023).

A later extension replaces the continuous family of jumps by a finite set of jump operators, possibly as few as one. The jump operators are built from user-chosen Hermitian couplings $Z=\mathrm{Tr}[e^{-\beta H}]$0, weighting functions $Z=\mathrm{Tr}[e^{-\beta H}]$1 supported on a finite interval, and inverse Fourier transforms $Z=\mathrm{Tr}[e^{-\beta H}]$2 with super-polynomial decay when $Z=\mathrm{Tr}[e^{-\beta H}]$3 is chosen from a compact-support Gevrey class. The resulting time-domain integrals can be truncated at $Z=\mathrm{Tr}[e^{-\beta H}]$4 and discretized efficiently (Ding et al., 2024).

3. Major algorithmic paradigms

One major family uses coherent block-encodings and purification. Chowdhury and Somma prepare the maximally entangled thermofield seed state, apply an approximation of $Z=\mathrm{Tr}[e^{-\beta H}]$5 using the Hubbard–Stratonovich transform and Hamiltonian simulation of $Z=\mathrm{Tr}[e^{-\beta H}]$6, and then use amplitude amplification. For local or sparse Hamiltonians this yields gate complexity

$Z=\mathrm{Tr}[e^{-\beta H}]$7

with polylogarithmic dependence on $Z=\mathrm{Tr}[e^{-\beta H}]$8 rather than polynomial dependence, and with a quadratic improvement in the dependence on $Z=\mathrm{Tr}[e^{-\beta H}]$9 over earlier approaches (1603.02940).

A distinct line aims at exact or near-exact sampling rather than approximate equilibration. Quantum “Coupling from the Past” adapts the Propp–Wilson perfect-sampling paradigm to energy eigenstates. It assumes an eigenbasis-preserving primitive quantum channel whose transition probabilities depend only on energies and uses phase estimation plus a classical coalescence structure to output perfect Gibbs samples without prior knowledge of the mixing time. The expected runtime depends strongly on degeneracy: for highly degenerate spectra it is polylogarithmic in the dimension and linear in the mixing time, while for non-degenerate spectra it is linear in the mixing time and quadratic in the dimension up to a logarithmic factor (França, 2017).

Metropolis-style constructions remain an active theme, but the technical route has shifted. A recent weak-measurement scheme replaces Marriott–Watrous rewinding by a weak accept/reject measurement and uses boosted quantum phase estimation based on medians of repeated QPE runs. One iteration defines a channel $H=\sum_i h_i$0 with expansion

$H=\sum_i h_i$1

where $H=\sum_i h_i$2 generates a continuous-time semigroup with a unique full-rank fixed point. If that semigroup mixes in time $H=\sum_i h_i$3, then $H=\sum_i h_i$4 discrete iterations suffice to reach $H=\sum_i h_i$5 up to controlled error. The construction was motivated in part by the observation that the earlier analysis of the Temme–Osborne–Vollbrecht–Poulin–Verstraete algorithm relied on a boosted and shift-invariant version of QPE that may not exist (Jiang et al., 2024).

The dissipative quantum Gibbs sampler of Zhang, Bosse, and Cubitt takes a different route. Its local step is the two-outcome instrument

$H=\sum_i h_i$6

with $H=\sum_i h_i$7 and a Hermitian step operator $H=\sum_i h_i$8 that factorizes into sequential weak measurements of the local Hamiltonian terms. The crucial point is that the Gibbs state is not obtained as the fixed point of a Markov process. Instead, the output is taken at a random stopping time $H=\sum_i h_i$9, and for a suitable stopping rule the expected output satisfies

$D$0

which approximates the Gibbs state. The paper states explicitly that, in contrast to the classical case, the quantum Gibbs state is not generated by converging to the fixed point of a Markov process, but by the states generated at the stopping time of a conditionally stopped process (Zhang et al., 2023).

4. Structure-exploiting constructions

Several samplers derive their efficiency from Hamiltonian structure rather than from a generic noncommutative template. For commuting local Hamiltonians, a recent reduction maps Gibbs-state preparation to classical Gibbs sampling. If $D$1 is a 2-local qudit commuting Hamiltonian, it can be mapped in $D$2 time to a 2-local classical Hamiltonian $D$3 on the same graph. If $D$4 is a 4-local qubit commuting Hamiltonian on a 2D lattice with no classical qubits, puncturing and a structure-lemma decomposition reduce the problem to a constant-local classical Hamiltonian on a planar graph. The corresponding quantum algorithm samples from $D$5, prepares the associated eigenstate of $D$6, and applies sequential randomized corrections; for the defected Toric code at any non-zero temperature, the total runtime is $D$7 (Hwang et al., 2024).

For smooth periodic potentials on the torus $D$8, a different paradigm leverages Fourier analyticity. The potential $D$9 is assumed $\rho$0 and analytically continuable to a complex strip, so its Fourier coefficients decay exponentially. Using zeroeth-order quantum queries, quantum Fourier transforms, and a quantum algorithm for linear ordinary differential equations, one solves a discretized Fokker–Planck equation whose stationary solution is the Gibbs density. The total gate cost is

$\rho$1

so the precision dependence is logarithmic, although the underlying Langevin mixing time can still be exponentially long in $\rho$2 for non-convex $\rho$3 (Motamedi et al., 2022).

Another alternative is quasi-probabilistic cluster expansion. Here the full Gibbs state is expanded into tensor products of local “Gibbs-cumulant” type states $\rho$4, which need not be positive semidefinite. Writing

$\rho$5

one samples cluster partitions according to positive weights, prepares the corresponding local mixed-state factors, and reweights measurements by an overall sign and a negativity factor $\rho$6. The method was demonstrated on 4-spin and 8-spin XY chains, including dynamical $\rho$7 correlations and specific heat (Eassa et al., 2023).

Adiabatic parent-Hamiltonian methods provide yet another viewpoint. Starting from a classical reversible Markov chain with detailed balance relative to a classical Gibbs distribution, one defines

$\rho$8

whose unique ground state encodes the Gibbs amplitudes. Speedups arise when the adiabatic path is detoured to avoid unfavorable critical points: in the 1D ferromagnetic Ising chain, detoured paths yield $\rho$9 rather than $\|\rho-\rho_\beta\|_1\le \epsilon$0; for weighted independent sets on certain graphs, one-parameter paths can scale as $\|\rho-\rho_\beta\|_1\le \epsilon$1 while detoured paths yield $\|\rho-\rho_\beta\|_1\le \epsilon$2. The weighted independent-set construction is naturally implementable on Rydberg atom arrays (Wild et al., 2020).

5. Complexity, mixing times, and lower bounds

Complexity analyses differ sharply by model class. For the dissipative stopped-process sampler, the expected output error obeys

$\|\rho-\rho_\beta\|_1\le \epsilon$3

while the stopping time satisfies an upper bound that leads, after parameter choice, to

$\|\rho-\rho_\beta\|_1\le \epsilon$4

times local-measurement cost $\|\rho-\rho_\beta\|_1\le \epsilon$5. The same work proves a fault-resilience statement: if each application of $\|\rho-\rho_\beta\|_1\le \epsilon$6 is implemented noisily with per-step error $\|\rho-\rho_\beta\|_1\le \epsilon$7, then the final sampling error is $\|\rho-\rho_\beta\|_1\le \epsilon$8 and does not grow with $\|\rho-\rho_\beta\|_1\le \epsilon$9 (Zhang et al., 2023).

For exact KMS-balanced noncommutative samplers, the implementation cost is organized around mixing time. The Chen–Kastoryano–Gilyén construction requires total Hamiltonian-simulation time

$$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$0

and a finite-jump-set variant yields total Hamiltonian-simulation time

$$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$1

with energy resolution depending only logarithmically on precision and mixing time (Chen et al., 2023, Ding et al., 2024).

Mixing-time upper bounds are only one side of the story. A generic bottleneck lemma for quantum Gibbs samplers proves exponential lower bounds whenever orthogonal high-mass sectors are separated by a low-mass bottleneck and the sampler has bounded Bohr-spectrum range or bounded locality range. For classical Hamiltonians with $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$2, the corresponding quantum Gibbs samplers also satisfy $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$3. The paper gives explicit families: random $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$4-SAT at low temperature with $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$5, $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$6-spin glasses with $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$7, good $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$8-qubit stabilizer codes with $$\rho_\beta=\frac{e^{-\beta H}}{\mathrm{Tr}[e^{-\beta H}]},$$9, and the ferromagnetic 2D transverse-field Ising model with $T=\beta^{-1}$0 for strictly local Lindblad operators (Gamarnik et al., 2024).

Recent work also studies acceleration relative to the Lindbladian gap $T=\beta^{-1}$1. For KMS-symmetric samplers with an explicit parent-Hamiltonian factorization $T=\beta^{-1}$2, purified Gibbs-state preparation can be reduced to a singular-value filtering problem. The resulting walk-free QSVT algorithm uses

$T=\beta^{-1}$3

queries, giving a quadratic improvement in spectral-gap dependence over $T=\beta^{-1}$4-type behavior for a broad class of samplers (Leng et al., 24 Apr 2026).

By contrast, there are settings where rigorous positive-gap guarantees are available. For truncated Coulomb gases and molecular systems in $T=\beta^{-1}$5, a quantum Markov semigroup tailored to the finite-rank truncation has a strictly positive spectral gap for every truncation, implying exponential convergence to the target Gibbs state. The paper presents this as the first rigorous mixing-time guarantee for Gibbs sampling in a Coulomb interacting continuous-variable quantum system (Becker et al., 16 Apr 2026).

6. Estimation tasks, hardware realizations, and computational complexity of sampling

Quantum Gibbs samplers are used not only for state preparation but also for thermodynamic estimation. In the dissipative stopped-process sampler, the stopping statistics of $T=\beta^{-1}$6 independent runs produce an unbiased estimator

$T=\beta^{-1}$7

with multiplicative $T=\beta^{-1}$8 error and variance $T=\beta^{-1}$9. For Coulomb gases and molecules, Gibbs sampling is embedded into a thermodynamic-integration procedure for estimating the free energy $\beta$00 with an end-to-end quantum algorithm whose complexity depends polynomially on $\beta$01, $\beta$02, and the inverse sampler gap (Zhang et al., 2023, Becker et al., 16 Apr 2026).

Experimental realizations emphasize a different regime. Studies of D-Wave quantum annealers report that programmable annealers behave as samplers generating independent configurations from low-temperature noisy Gibbs distributions. On up to 16 qubits with native Chimera connectivity and $\beta$03, a “sweet-spot”

$\beta$04

minimizes both control-noise and residual-transverse-field distortion; across 13 models, the total-variation distance between the empirical distribution and the best-fit Gibbs law drops below $\beta$05 on every instance and often reaches $\beta$06–$\beta$07 (Vuffray et al., 2020, Nelson et al., 2021).

The computational-complexity side has developed in parallel. For carefully chosen commuting parent Hamiltonians of shallow circuits, thermalization under Davies-type Lindbladians mixes rapidly at constant temperature, yet sampling from the measurement distribution of the Gibbs state is classically hard under standard complexity assumptions. One construction proves quantum computational advantage for commuting local Hamiltonians at constant temperature, with mixing time $\beta$08 and quantum runtime $\beta$09 when $\beta$10 (Bergamaschi et al., 2024). A subsequent strengthening shows that hardness persists even for 5-local Hamiltonians on a 3D cubic lattice and for 6-local Hamiltonians with constant additive sampling error, while polynomial-time quantum preparation remains possible; the hardness is also robust to imperfect measurements (Rajakumar et al., 2024).

A recurrent misunderstanding is that finite temperature automatically renders sampling classically benign. The constant-temperature hardness results show otherwise for specific families, while the slow-mixing results show that locality and detailed balance alone do not guarantee efficient equilibration. Conversely, the commuting-Hamiltonian reductions, continuous-potential algorithms, and hardware studies show that substantial positive results do exist, but only under structural assumptions that must be stated explicitly.