Direct Estimation of the Density of States for Fermionic Systems

TL;DR

This paper introduces quantum algorithms to directly estimate the density of states (DOS) and derived thermodynamic properties for quantum systems, with a focus on fermionic subspaces. It combines subspace DOS estimation, random-state initialization, and Gaussian-windowed Fourier reconstruction to enable DOS recovery with shallow, noisy circuits, while remaining compatible with variational approaches for NISQ devices. The authors demonstrate robustness against algorithmic and gate noise, provide detailed numerical experiments on Fermi-Hubbard and spin models, and discuss resource implications and extensions toward fault-tolerant quantum computation. Overall, the work offers a practical pathway toward early quantum advantage in computing thermodynamic quantities for complex quantum materials and chemistry problems.

Abstract

Simulating time evolution is one of the most natural applications of quantum computers and is thus one of the most promising prospects for achieving practical quantum advantage. Here, we develop quantum algorithms to extract thermodynamic properties by estimating the density of states (DOS), a central object in quantum statistical mechanics. We introduce several innovations that significantly improve the practicality and extend the generality of previous techniques. First, our approach allows one to estimate the DOS only for a specific subspace of the full Hilbert space. This is crucial for fermionic systems, since both canonical and grand canonical ensemble thermal equilibrium properties depend on subspaces of fixed number. Second, in our approach, by time evolving very simple, random initial states, such as randomly chosen computational basis states, we can exactly recover the DOS on average. Third, due to circuit-depth limitations, we only reconstruct the DOS up to a convolution with a Gaussian window - thus all imperfections that shift the energy levels by less than the width of the convolution window will not significantly affect the estimated DOS. For these reasons we find the approach is a promising candidate for early quantum advantage as even short-time, noisy dynamics can yield a semi-quantitative reconstruction of the DOS (convolution with a broad Gaussian window), while early fault tolerant devices will likely enable higher resolution DOS reconstruction through longer time evolutions. We demonstrate the practicality of our approach in representative Fermi-Hubbard and spin models and indeed find that our approach is highly robust against algorithmic errors in the time evolution and against gate noise. We further demonstrate that our approach is compatible with NISQ-friendly variational techniques, introducing and leveraging a new technique for variational time evolution.

Figures (12)

• Figure 1: Density of states of a Hubbard model on subspaces of fixed number $M$. We cumulatively plot $g_M(E)$ (up to a broadening of lines discussed in \ref{['sec:window_width_effects']}) for a $(3\times 2)$ open-boundary grid Hubbard model \ref{['eqn:fermi_hubbard_hamiltonian']} with $J=-1$, $U=2$ for a series of different particle numbers $M$. Each $M$ has a different DOS (different colours) and a different normalization factor $|\mathcal{S}|$ which has important implications for computing the partition function $Z_M(\beta)$ or the grand canonical partition function $\mathcal{Z}(\beta,\mu)$. Here, $M=6$ is the most populous subspace.
• Figure 2: Quantum circuits for evaluating the FDOS.(a) DQC1 circuits: First, setting $\mathcal{V}$ to be the Hadamard transform $H^{\otimes n}$ allows us to estimate the FDOS over the entire Hilbert space as relevant for \ref{['sec:constant_particle_numb']}. This way Bell pairs are created between the two registers such that the reduced density matrix of the lower, main register is maximally mixed -- as in standard DQC1 computations. Second, setting $\mathcal{V}$ to be a circuit that prepares a uniform superposition $|\mathcal{S}|^{-1/2} \sum_k \ket{\psi_k}$ allows us to estimate the FDOS $G_\mathcal{S}(t)$ on a subspace (\ref{['sec:variable_particle_numb']}) spanned by $\{\ket{\psi_k}\}_{k\in \mathcal{S}}$. We detail an efficient construction of $\mathcal{V}$ for fermionic Hamiltonians in \ref{['sec:dqc_subspace']}. (b) Alternatively, the FDOS can be estimated by random state sampling methods outlined in \ref{['sec:dos_with_random_states']}. Here the main register is initialised in a randomly chosen initial state and and thus the need for a second register to prepare a maximally mixed state is eliminated entirely. The random intialisation circuit $\Psi$ can be extremely simple - even random single-qubit bit-flips suffice, but we explore further variants in \ref{['sec:state_sampling_numerics']}. The real (imaginary) parts of the FDOS are encoded as the ancilla probability by setting $W{=}H$ ($W{=}S^\dagger H$).
• Figure 3: Effects of finite simulation length. We study a $(3\times2)$ grid Fermi-Hubbard model with open boundary conditions, nearest-neighbour couplings, $J=-1$ and $U=2$, and assuming the fixed-number subspace $M=6$. (a) The FDOS $G_6(t)$ (black) is estimated using a quantum computer and allows us to reconstruct the DOS $g_6(E)$ on a subspace of fixed particle number $M=6$ through inverse Fourier transforming. As detailed in \ref{['sec:window_width_effects']}, to cap the total time evolution length, the signal is multiplied in post-processing by a Gaussian Fourier window $W_\sigma(t)$ of temporal width $\sigma=180$ (blue), $\sigma=60$ (orange), and $\sigma=20$ (green) -- shaded regions represent the amplitude of the window function. (b) different resolution-limited DOS $\tilde{g}^{\sigma}(E)$ are obtained from each signal in (a): wider temporal window width (blue) leads to finer-resolution spectral features while shorter time evolutions (green) blur out spectral features. (c) the partition function $Z_6(T)$ from \ref{['eqn:canonical_partition_function']} is estimated for increasing temperatures $T=1/\beta$ using the resolution-limited DOS $\tilde{g}^{\sigma}(E)$ from (b) -- the shaded grey region represents the true partition function. Low-temperature physics imposes stricter requirements on the duration of dynamics: shorter temporal window widths (green and orange) yield poorer approximations at low temperatures due to the concentration of the Boltzmann factor $e^{-\beta E}$ at low energies. In contrast, the inset illustrates that even short-time dynamics can be sufficient for accurate prediction of higher temperature properties.
• Figure 4: Comparison of random initial-state sampling methods. Error $\epsilon(f,\sigma)$ as a function of per-timestep shot budget $N_s$ in the DOS of a Heisenberg model in \ref{['eqn:heisenberg_hamiltonian']} using different random sampling methods (Haar-random, single-qubit continuous rotations, and random bit flips). As in most hardware platforms loading a new circuit is expensive and thus multiple repetitions of a single circuit is desirable, (coloured lines) represent that the same random initial state was reused $N_r$ times. Lower dark grey shaded regions represent the errors obtained by DQC1 computation via \ref{['fig:circ']}(a) whereas the light grey shaded area represents completely random noise. Parallel lines in the above log-log plots confirm standard shot-noise scaling $N_s = O(\epsilon^{-2})$ from \ref{['stat:shot_number']} for all techniques, i.e., performances are only different by constant absolute factors, and for $N_r=1$ the random-sampling approach has the same sample complexity as the DQC1 computation as expected. For larger values of $N_r$ (repeated initial samples), the more uniform random sampling methods perform slightly better, however, they may require more complex initialisation circuits -- the single-qubit, random Euler angle approach is likely to be the most practical one. Confidence intervals are too small to visualize on this plot.
• Figure 5: Effect of initial-state-preparation circuit depth on DOS error.(a) We consider here random initial state sampling in \ref{['fig:circ']} through randomly choosing gate parameters in the pictured layered ansatz. Each layer consists of arbitrary single-qubit rotations (parametrized by Tait-Bryan angles) and nearest-neighbour $ZZ$ rotations. (b) Error $\epsilon(f,\sigma)$ in the estimated DOS for an increasing number of layers $L$, at three different sampling circuit repetitions $N_r\in\{1,10,10^2\}$. Shaded bars represent bootstrapped 95% confidence intervals. Increasing $L$ interpolates between the $1$-design and $2$-design limits, the latter of which leads to reduced impact of sampling circuit repetition. Indeed more uniform random sampling of the Hilbert space is preferable when reusing initialisation circuits (via $N_r >1$).

Theorems & Definitions (3)

• proof
• proof
• proof