XDiag: Exact Diagonalization for Quantum Many-Body Systems

TL;DR

XDiag presents an open-source exact diagonalization toolkit for quantum many-body systems that fuses symmetry-adapted bases, sublattice coding, Lin tables, and random-hashing with a high-performance C++ core and a Julia wrapper. The work introduces the first public implementation of sublattice coding for large-scale spin diagonalizations and demonstrates near-linear scaling on thousands of CPU cores across shared- and distributed-memory configurations, while supporting multiple Hilbert space types such as $S=1/2$ spins, Hubbard, and $t$-$J$ models. It provides extensive documentation, a user guide, 20+ examples (ground-state to thermal states), and reproducible benchmarks, highlighting the toolkit’s versatility for ground-state, spectral, dynamical, and thermodynamic analyses. The combination of optimized algorithms, advanced symmetry handling, and a user-friendly scripting interface offers a powerful platform for high-precision quantum many-body simulations with broad practical impact in computational physics.

Abstract

Exact diagonalization (ED) is a cornerstone technique in quantum many-body physics, enabling precise solutions to the Schrödinger equation for interacting quantum systems. Despite its utility in studying ground states, excited states, and dynamical behaviors, the exponential growth of the Hilbert space with system size presents significant computational challenges. We introduce XDiag, an open-source software package designed to combine advanced and efficient algorithms for ED with and without symmetry-adapted bases with user-friendly interfaces. Implemented in C++ for computational efficiency and wrapped in Julia for ease of use, XDiag provides a comprehensive toolkit for ED calculations. Key features of XDiag include the first publicly accessible implementation of sublattice coding algorithms for large-scale spin system diagonalizations, efficient Lin table algorithms for symmetry lookups, and random-hashing techniques for distributed memory parallelization. The library supports various Hilbert space types (e.g., spin-1/2, electron, and t-J models), facilitates symmetry-adapted block calculations, and automates symmetry considerations. The package is complemented by extensive documentation, a user guide, reproducible benchmarks demonstrating near-linear scaling on thousands of CPU cores, and over 20 examples covering ground-state calculations, spectral functions, time evolution, and thermal states. By integrating high-performance computing with accessible scripting capabilities, XDiag allows researchers to perform state-of-the-art ED simulations and explore quantum many-body phenomena with unprecedented flexibility and efficiency.

Figures (6)

• Figure 1: Enumeration of computational basis states of different Hilbert space blocks. (a) block on $N=4$ sites with $n_{\uparrow}=2$$\uparrow$-spins. (b) block on $N=4$ sites with $n_{\uparrow}=2$$\uparrow$-spins and $n_{\downarrow}=1$$\downarrow$-spins. (c) block on $N=4$ sites with $n_{\uparrow}=2$$\uparrow$-spins and $n_{\downarrow}=1$$\downarrow$-spins. The index of a given basis configuration (represented as an object of type ) can be obtained using the function.
• Figure 2: Enumeration of computational basis states of symmetry-adapted Hilbert space blocks with translational symmetry in the $k=0$ representation. (a) block on $N=4$ sites with $n_{\uparrow}=2$$\uparrow$-spins. (b) block on $N=4$ sites with $n_{\uparrow}=2$$\uparrow$-spins and $n_{\downarrow}=1$$\downarrow$-spins. (c) block on $N=4$ sites with $n_{\uparrow}=2$$\uparrow$-spins and $n_{\downarrow}=1$$\downarrow$-spins.
• Figure 3: Scaling of the computation time for a single Lanczos iteration as a function of the total number of threads used. The calculations were done for the Heisenberg $a)$, Hubbard $b)$ and t-J $c)$ models on a linear chain using different total numbers of spins/sites. The dashed black line in each panel represents the slope corresponding to ideal linear scaling. In the calculation, we used both the U(1) and lattice symmetries.
• Figure 4: Scaling of the computation time for a single Lanczos iteration as a function of the total number of MPI processes used. The calculations were done for the Heisenberg $a)$, Hubbard $b)$, and t-J $c)$ models using different total numbers of spins/sites. The dashed black line in each panel represents the slope corresponding to ideal linear scaling. In $a)$, the calculation was performed in the zero magnetization sector, in $b)$ it was done at half-filling for $N =20$, at $3/11$ for $N=22$ and $1/3$ for $N=24$.
• Figure 5: Computational time for the generation of the Hamiltonian matrix in the sparse CSR matrix format as a function of the total number of threads used. The calculations were done for the Heisenberg $a)$, Hubbard $b)$, and t-J $c)$ models using different total numbers of spins/sites. The dashed black line in each panel represents the slope corresponding to ideal linear scaling.