We present QCommute, a software tool implemented in C++ for symbolic computation of nested commutators between a Hamiltonian and local observables in quantum many-body spin-1/2 systems on one-, two-, and three-dimensional hypercubic lattices. The computation is performed algebraically directly in the thermodynamic limit, and the Hamiltonian parameters are kept symbolic. Importantly, this way the entire parameter space is covered in a single run. The implementation supports extensive parallelization to achieve high computational performance. QCommute enables the investigation of quantum dynamics in strongly correlated regimes that are inaccessible to perturbative approaches, either through direct Taylor expansion in time or via advanced techniques such as the recursion method.
We present QCommute, a software tool implemented in C++ for symbolic computation of nested commutators between a Hamiltonian and local observables in quantum many-body spin- 1/2 1 / 2 systems on one-, two-, and three-dimensional hypercubic lattices. The computation is performed algebraically directly in the thermodynamic limit, and the Hamiltonian parameters are kept symbolic. Importantly, this way the entire parameter space is covered in a single run. The implementation supports extensive parallelization to achieve high computational performance. QCommute can serve as a computational backend for Heisenberg-picture approaches to quantum dynamics in strongly correlated regimes, ranging from direct Taylor expansion in time to advanced techniques such as the recursion method.
We report an implementation of the recursion method that addresses quantum many-body dynamics in the nonperturbative regime. The method essentially amounts to constructing a Lanczos basis in the space of operators and solving coupled Heisenberg equations in this basis. The reported implementation has two key ingredients: a computer-algebraic routine for the symbolic calculation of nested commutators and a procedure to extrapolate the sequence of Lanczos coefficients according to the universal operator growth hypothesis. We apply the method to calculate infinite-temperature correlation functions for spin-
1/2 systems on one- and two-dimensional lattices. In two dimensions the accessible timescale is large enough to essentially embrace the relaxation to equilibrium. The method allows one to accurately calculate transport coefficients. As an illustration, we compute the diffusion constant for the transverse-field Ising model on a square lattice.
High-performance numerical methods are essential not only for advancing quantum many-body physics but also for enabling integration with emerging quantum computing platforms. We present a scalable, general-purpose parallel algorithm for quantum simulations based on -represented quantum algebra (ORQA). This framework applies to arbitrary spin systems and naturally integrates with quantum circuit simulation in the Heisenberg picture, making it particularly relevant to recent large-scale experiments on superconducting qubit processors [Y. Kim ]. As a benchmark, we simulate the kicked Ising model on a 127-qubit heavy-hexagon lattice, tracking the time evolution of local magnetization using up to 1 × 10 12 Pauli strings. Executed on the supercomputer Fugaku, our simulations exhibit strong scaling up to 2 17 parallel processes with near-linear communication overhead. These results establish ORQA as a practical and high-performance tool for quantum many-body dynamics and highlight its potential for integration into hybrid quantum-classical computational frameworks, complementing recent advances in tensor network and surrogate simulation techniques.
We present a symbolic implementation of the recursion method for dynamical correlations and transport in fermionic systems on one-, two-, and three-dimensional lattices. The implementation is applicable in the strongly correlated regime, yields results directly in the thermodynamic limit, and covers all time scales, from short times through the thermalization time to the late-time asymptotics. Focusing on two paradigmatic models -- interacting spinless fermions and the Hubbard model -- we confirm the universal operator growth hypothesis in the fermionic case, compute infinite-temperature current-current autocorrelation functions, and determine the charge diffusion constant and the high-temperature conductivity. The diffusion constant is obtained in both the perturbative and the nonperturbative regime of the interaction strength, and we locate the boundary between them. We compare our approach with the Majorana propagation method, which converges up to intermediate times and agrees with ours where converged. Our results highlight a symbolic computational paradigm in which the most expensive step is performed once, producing reusable symbolic output that yields physical observables for arbitrary model parameters.
The recursion method, which solves coupled Heisenberg equations in a Lanczos operator basis, has recently emerged as a powerful nonperturbative tool for computing dynamical correlation functions in strongly correlated two- and three-dimensional quantum many-body systems. Motivated by this success, we investigate whether the method can be extended to expectation values of observables following a quantum quench. We find that such an extension encounters an obstacle absent in the computation of dynamical correlation functions. The latter are fully determined by the Lanczos coefficients b_n b n , which in generic systems exhibit universal behavior, enabling reliable extrapolation from the first few dozens of explicitly computed coefficients. In contrast, quench dynamics additionally requires “quench coefficients” c_n c n , defined as overlaps of Lanczos basis operators with the initial state. We show that, unlike the Lanczos coefficients, the quench coefficients exhibit no universal structure and cannot be reliably extrapolated, thereby limiting the time up to which the method yields accurate results. The behavior of quench coefficients is highly state-dependent, ranging from decaying to irregular or even growing sequences; typically, the less regular the sequence c_n c n , the shorter the accessible timescale. Nevertheless, for favorable initial states, the method remains competitive with state-of-the-art approaches. Moreover, its symbolic implementation allows a single computation to be reused across different Hamiltonian parameters and initial states, making it particularly advantageous in studies requiring extensive scans over Hamiltonian parameters or initial states.
This thesis focuses on the development, application, and high-performance implementation of numerical methods for simulating quantum many-body systems. In the first chapter, I provide an overview of several most successful conventional approaches for simulating quantum many-body systems, along with their current limitations. Then, the basic properties of the quantum spin model, which will be the subject of study in the following chapters, and its ground state are briefly summarized. Additionally, I discuss the background and the motivations of the research projects carried out during my PhD research. Chapter 2 introduces two numerical algorithms that are extensively used in this thesis, which are exact diagonalization (ED) and variational Monte Carlo. I outline the process of constructing a Hamiltonian matrix and employing iterative methods for exact diagonalization. Subsequently, I introduce the mechanism of the four important components in the ``standard" variational Monte Carlo that are widely used in the most recent research. In Chapter 3, we propose a symmetry-projected variational ansatz in the form of linear combinations of simple restricted Boltzmann machines. This construction allows one to explore states outside of the original variational manifold and increase the representation power with moderate computational effort. Besides allowing one to restore spatial symmetries, an expansion in terms of Krylov states using a Lanczos recursion offers a solution that can further improve the quantum state accuracy. We illustrate these ideas with an application to the Heisenberg
J1−J2 model on the square lattice, a paradigmatic problem under debate in condensed matter physics,and achieve state-of-the-art accuracy in the representation of the ground state. In Chapter 4, we present a variational approach based on a Chebyshev expansion of the spectral function, and each Chebyshev polynomial is represented by a neural network quantum state. The Chebyshev moments are obtained by recursively applying the Hamiltonian and projecting on the space of variational states using a modified natural gradient descent method. We compare this approach with a modified approximation of the spectral function which uses a Krylov subspace constructed from the "Chebyshev wave functions." We present results for the one-dimensional and two-dimensional Heisenberg models on the square lattice and compare them with those obtained by other methods in the literature. In Chapter 5, we describe a method to carry out thermodynamic simulations of quantum many-body models using typicality and variational representations of quantum states. Minimally entangled typical thermal states (METTS) are a construction that allows one to solve for the imaginary time evolution of quantum many-body systems. By using wave functions that are weakly entangled, one can take advantage of efficient representations in the form of matrixproduct states. We generalize these ideas to arbitrary variational wave functions and we focus, as an illustration, on the particular case of restricted Boltzmann machines. We demonstrate the application of the method to quantum spin systems in one and two spatial dimensions. In Chapter 6, we study the properties of the impurity band in heavily doped nonmagnetic semiconductors using the Jacobi-Davidson algorithm and the supervised deep learning method. The disorder averaged inverse participation ratio and Thouless number calculation show us the rich structure inside the impurity band. A convolutional neural network model, which is trained to distinguish the extended/localized phase of the Anderson model with high accuracy, shows us the results in good agreement with the conventional approach. Together, we find that there are three mobility edges in the impurity band for a specific on-site impurity potential, which means the presence of the extended states while filling the impurity band. In the final chapter, Chapter 7, we discuss the implementation of architecture-aware optimizations for spin dynamics simulations on x86 CPUs and Nvidia GPUs. The time cost of evaluating the spin-spin correlation functions in real space increases quadratically as the number of spins N, leading to significant computational effort, making the simulation of large spin systems very challenging. We propose to use a highly optimized general matrix multiply (GEMM) subroutine to calculate the dynamical spin-spin correlation function that can achieve near-optimal hardware utilization. Furthermore, we fuse the element-wise operations in the calculation of S(q,t) into the in-house GEMM kernel to reduce memory latency. On Nvidia GPU, we perform the on-the-fly calculation in the epilogue of the GEMM subroutine to avoid saving intermediate results on global memory, which significantly reduces memory usage and achieves better performance than
cuBLAS.--Author's abstract