A symmetry-based guarantee for simulating interacting fermions

What the researchers establish

In an arXiv preprint, Chaithanya Rayudu, Jun Takahashi, Ojas Parekh and Andrew Zhao present a classical randomized algorithm for estimating the partition function of a specified family of interacting fermion systems. Fermions include electrons, and Hubbard models provide simplified descriptions of how their motion and interactions combine.

The partition function sums the statistical weights of a system’s possible energy states at a given temperature. It determines the free energy, and sufficiently low-temperature estimates can also approximate the ground energy—the lowest energy the system can have. Computing it directly becomes difficult as the number of fermionic modes grows.

The authors prove that their algorithm estimates this quantity to a chosen relative error in polynomial time. The bound depends on the number of modes, the number of interaction terms, inverse temperature, an energy-scale parameter and inverse target error. Their stated success probability is at least three quarters. They also derive polynomial-time guarantees for estimating free energy and ground energy to a chosen additive accuracy.

Which systems qualify?

The result applies to what the paper calls the attractive Kramers class. Its interactions attract fermions occupying paired modes related by a particular time-reversal symmetry. Mathematically, the relevant time-reversal operation squares to minus the identity, and the single-particle motion respects that symmetry.

Examples include attractive Fermi-Hubbard models with real hopping amplitudes on any graph. The result also covers repulsive Fermi-Hubbard models at half filling, with real hopping amplitudes on bipartite graphs, through a staggered particle-hole transformation. Bipartite means the sites can be divided into two groups with connections only between groups. These qualifications matter: the theorem does not cover repulsive Hubbard models on bipartite graphs generally.

The guarantee does not require weak interactions. However, runtime still depends polynomially on the specified energy scale, which includes the interaction strengths. Strong-interaction coverage therefore does not mean that increasing those strengths has no computational cost.

How the algorithm works

Quantum Monte Carlo methods turn a quantum calculation into a statistical sampling problem. For fermions, individual contributions can be negative or complex. Cancellations then make accurate estimation difficult—the notorious sign problem. Even when all contributions are nonnegative, a sampler can still take too long to reach its target distribution.

The researchers first divide inverse temperature into many small steps, using a controlled approximation called Trotterization. They expand the interaction factors into classical configurations. Wick’s theorem, a standard identity for noninteracting fermions, expresses each configuration’s weight using a determinant of a matrix built from two-point correlations.

Kramers symmetry gives that matrix a special structure, called quaternionic. The relevant determinants are consequently nonnegative. This supplies usable sampling weights without the sign problem.

The next step is the distinctive one: a determinantal worm Markov chain. Rather than sampling only the configurations needed for the final answer, the algorithm temporarily explores a larger space containing auxiliary defects. In physical terms, the enlargement corresponds to introducing a weak field that creates or annihilates fermion pairs. Configurations without those defects recover the desired distribution.

The authors prove that a suitable field keeps defect-free configurations sufficiently common: their probability is at least one half. Thus, useful samples do not become prohibitively rare.

Why the sampler converges—and what remains unproven

The convergence proof uses strong log-concavity, a mathematical property of a polynomial encoding the configuration weights. Combined with an existing result about reweighted sampling distributions, this yields a polynomial bound on the sampler’s mixing time: how long it takes to approach its target distribution. Sampling at a sequence of interaction strengths then provides the partition-function estimate.

This is a theoretical algorithm and proof, not a reported benchmark showing faster calculations on particular materials. Polynomial scaling alone does not establish practical competitiveness. Nor does the analysis prove fast convergence for every existing determinant Monte Carlo implementation; the authors’ sampler is different.

The covered systems also possess a Lee-Yang property, restricting where partition-function zeros can occur when an external field is allowed to be complex. The authors suggest that combining such restrictions with freedom from the sign problem could guide broader efficient algorithms. That is a research direction, not a theorem for every system sharing those properties. The paper explicitly distinguishes its covered class from broader sign-problem-free Kramers formulations.