Scalable Quantum Monte Carlo Method for Polariton Chemistry via Mixed Block Sparsity and Tensor Hypercontraction Method

We present a reduced-scaling auxiliary-field quantum Monte Carlo (AFQMC) framework designed for large molecular systems and ensembles, with or without coupling to optical cavities. Our approach leverages the natural block sparsity of Cholesky decomposition (CD) of electron repulsion integrals in molecular ensembles and employs tensor hypercontraction (THC) to efficiently compress low-rank Cholesky blocks. By representing the Cholesky vectors in a mixed format, keeping high-rank blocks in block-sparse form and compressing low-rank blocks with THC, we reduce the scaling of exchange-energy evaluation from quartic to robust cubic in the number of molecular orbitals, while lowering memory from cubic toward quadratic. Benchmark analyses on one-, two-, and three-dimensional molecular ensembles (up to ~1,200 orbitals) show that: a) the number of nonzeros in Cholesky tensors grows linearly with system size across dimensions; b) the average numerical rank increases sublinearly and does not saturate at these sizes; and (c) rank heterogeneity-some blocks nearly full rank and many low rank, naturally motivating the proposed mixed block sparsity and THC scheme for efficient calculation of exchange energy. We demonstrate that the mixed scheme yields cubic CPU-time scaling with favorable prefactors and preserves AFQMC accuracy.

💡 Research Summary

**
The manuscript introduces a reduced‑scaling auxiliary‑field quantum Monte Carlo (AFQMC) framework tailored for large molecular ensembles, both in the presence and absence of optical cavity coupling. The central computational bottleneck in AFQMC is the evaluation of the two‑electron exchange term, which traditionally scales as O(N⁴) with the number of molecular orbitals N. The authors address this by exploiting two structural properties of the Cholesky‑decomposed electron‑repulsion integrals (ERIs) that arise in molecular ensembles: (i) block sparsity due to spatial locality and intermolecular separation, and (ii) low‑rank character of many Cholesky blocks.

First, the Pauli‑Fierz Hamiltonian is recast into the Monte‑Carlo Hamiltonian form required by AFQMC. The ERIs are expressed via a Cholesky decomposition V_{pqrs}=∑_γ L_γ(pq)L_γ(rs), where the number of Cholesky vectors N_γ grows linearly with N. In extended ensembles, the Cholesky tensors naturally form a block‑sparse pattern: in one dimension they are tridiagonal, while in two and three dimensions each block connects only to a constant number of neighboring blocks (degree d). Consequently, the total number of non‑zero elements in each L_γ scales linearly with N, and the storage requirement for the half‑rotated tensors is O(O·N·N_γ) (O denotes the number of occupied orbitals).

Second, the authors perform a detailed rank analysis of each block. While a minority of blocks (typically those with the largest norms) retain a near‑full rank close to N, the majority exhibit substantially lower numerical rank R_γ. For low‑rank blocks, tensor hypercontraction (THC) provides an efficient factorization L_γ(pq)≈∑_μ X_γ(p,μ)U_γ(q,μ), reducing storage from O(N²) to O(N·R_γ) and the cost of forming the intermediate f_γ^{ij}=∑_p L_γ(p,i)Θ(p,j) from O(N²) to O(N·R_γ).

The key innovation is a mixed representation that stores high‑rank blocks in their original block‑sparse form and compresses low‑rank blocks with THC. An analytical decision rule is derived by equating the per‑γ cost of the two approaches, yielding a size‑independent rank threshold R*≈κ(d+1)s (s is the block size, κ absorbs constant prefactors). Blocks with R_γ≤R* are treated with THC, while the rest remain block‑sparse. This hybrid scheme leads to an overall exchange‑energy evaluation cost that scales as ∼O(N³) with a memory footprint that scales roughly as O(N²), a substantial improvement over pure CD (O(N³) memory, O(N⁴) time) or pure THC (sub‑quartic time but large prefactors for realistic system sizes).

Benchmark calculations on one‑, two‑, and three‑dimensional molecular ensembles containing up to ~1,200 orbitals confirm the theoretical scaling. The number of non‑zero Cholesky elements grows linearly with system size, and the average numerical rank increases sub‑linearly without saturation at the examined sizes. The mixed BS‑THC approach delivers cubic wall‑time scaling with favorable prefactors and retains the high accuracy of standard AFQMC; energy differences relative to conventional AFQMC are within 10⁻³ Hartree for test polaritonic systems.

Furthermore, the electron‑photon coupling terms, which involve only a few photon modes (typically <10), add only O(N·P) cost, where P is the number of photon modes. Thus the method is especially advantageous for collective strong‑coupling regimes where many molecules interact coherently with a small set of cavity modes.

In summary, the paper presents a practical, scalable AFQMC algorithm that combines block sparsity and tensor hypercontraction to overcome the O(N⁴) exchange bottleneck. It enables accurate quantum‑chemical simulations of large polaritonic ensembles, opening the door to first‑principles studies of cavity‑modified chemistry and strongly correlated light‑matter systems that were previously out of reach.