SmoQyElPhQMC.jl: An open-source Julia package for efficient and scalable quantum Monte Carlo simulations of electron-phonon coupled models

Abstract: We introduce version 1.0 of the SmoQyElPhQMC package, an open-source Julia code for performing scalable quantum Monte Carlo simulations of electron-phonon coupled model Hamiltonians. SmoQyElPhQMC is built upon the SmoQyDQMC codebase and implements improved versions of the algorithms presented in [B. Cohen-Stead \textit{et al}., Phys. Rev. E {\bf 105}, 065302 (2022)] to enable linear-scaling simulations of a broad class of uncorrelated $e$-ph models both in system size and inverse temperature. By extending the functionality of the flexible scripting interface introduced in SmoQyDQMC, the SmoQyElPhQMC package continues to allow users to adapt it to different workflows and interface with other software packages in the Julia ecosystem. The code for this package can be downloaded from our GitHub repository at https://github.com/SmoQySuite/SmoQyElPhQMC.jl or installed using the Julia package manager. The online documentation, including examples, can be obtained from our documentation page at https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/.

• The paper introduces SmoQyElPhQMC.jl which implements near-linear scaling QMC techniques using hybrid Monte Carlo and exact Fourier acceleration for electron–phonon systems.
• It demonstrates significant performance improvements over traditional DQMC, accurately capturing order parameters and spectral features in benchmark Holstein and oSSH models.
• Its innovative preconditioning, block matrix approaches, and stochastic measurement methods enable precise and scalable studies of complex electron–phonon interactions.

SmoQyElPhQMC.jl: Efficient and Scalable Quantum Monte Carlo for Electron-Phonon Models

Introduction and Context

The SmoQyElPhQMC.jl package represents a significant advancement in large-scale quantum Monte Carlo (QMC) simulations of electron-phonon (e-ph) coupled models. Built atop the SmoQyDQMC.jl framework, it implements a suite of algorithmic optimizations that enable near-linear scaling in both system size $N$ and inverse temperature $\beta$, targeting a broad class of uncorrelated e-ph Hamiltonians. The software leverages hybrid Monte Carlo (HMC) strategies, exact Fourier acceleration (EFA), advanced preconditioning, and stochastic measurement techniques for efficient simulation, incorporating flexible support for generic tight-binding models, various phonon branches, long-range and nonlinear e-ph couplings, disorder, and dynamic carrier density control.

This development addresses a prominent computational bottleneck in existing QMC schemes for e-ph systems. Traditional DQMC-based approaches, while robust and widely adopted, exhibit cubic or worse scaling with system size—rendering them prohibitive in the adiabatic regime or for large $N$ and $\beta$ relevant to real materials. SmoQyElPhQMC.jl instead enables operability at unprecedented scale, accuracy, and model generality.

Supported Models and Hamiltonian Structure

The code supports an extensive variety of e-ph coupled Hamiltonians, formulated as

$$\hat{\mathcal{H}} = \hat{\mathcal{U}} + \hat{\mathcal{V}} + \hat{\mathcal{K}}$$

where $\hat{\mathcal{U}}$ describes phononic (harmonic, anharmonic, and dispersive) terms; $\hat{\mathcal{V}}$ includes the electron on-site potential and diagonal e-ph couplings (e.g., Holstein/Fröhlich types), and $\hat{\mathcal{K}}$ contains the electronic hopping (including phonon-modulated SSH/Peierls-like terms).

The package enables:

• Arbitrary lattice geometry and basis dimension,
• Momentum-dependent, long-range, and nonlinear e-ph interactions,
• Multiple phonon branches with independent microscopic coupling,
• Low-energy (optical/acoustic) phonons and anharmonic lattice potentials,
• Chemical potential feedback for fixed carrier density,
• Full spatial disorder support,
• Specialized observable measurement infrastructure.

Algorithmic Innovations

At its core, SmoQyElPhQMC.jl utilizes HMC—specifically, the EFA-HQMC method—for global phonon updates, effectively sampling the high-dimensional continuous phonon field with low autocorrelation and enabling ergodic exploration. The key algorithmic contributions are as follows:

• Checkerboard Exponential (MSCHK) Approximation: Sparse exponentiation of the kinetic term preserves scalability for arbitrary hopping range and generic SSH couplings while introducing only $\mathcal{O}(\Delta\tau^2)$ systematic error.
• Fermion Determinant as Pseudofermion Integral: The costliest operation in DQMC (evaluating $|\det M|^2$) is recast as a complex Gaussian integral over pseudofermion fields, transforming the bottleneck into an iterative linear solve amenable to preconditioning and parallelization.
• Block Structure and Preconditioning: The package replaces the conventional small-matrix representation with a block-sparse ($\beta$0) matrix for efficient storage and matrix-vector operations, using a preconditioner inspired by the adiabatic limit and implemented via the Kernel Polynomial Method (KPM) for optimal conditioning of the conjugate gradient (CG) solver.
• Ergodicity Enhancements: Reflection and swap updates are implemented to overcome topological barriers and ergodicity breaking, crucial for reliable sampling especially in models with discrete symmetry breaking.
• Stochastic Measurements: All Green's function and correlation observables are computed via stochastic estimators with $\beta$1 random (phase) vectors, ensuring unbiased results at $\beta$2 cost, with autocorrelation further reduced by FFT-based translational averaging.

These features result in a robust toolkit for nonperturbative, sign-problem-free simulation of large e-ph Hamiltonians across the full range of physically relevant parameters.

Numerical Benchmarks and Physical Results

Benchmark calculations focus on the square lattice Holstein and optical SSH (oSSH) models at half-filling in the adiabatic regime ($\beta$3), chosen for their direct relevance to charge/bond-order instabilities in low-dimensional electron-lattice systems.

Simulations quantitatively capture:

• The onset and scaling of charge density wave (CDW) and bond order wave (BOW) order parameters,
• Near-linear runtime scaling in both system size and inverse temperature, with only a weak residual dependence due to CG iteration counts,
• High acceptance rates ($\beta$4) for HMC trajectories,
• Accurate spectral functions $\beta$5, with robust signature of long-range order and spectral weight redistribution due to e-ph coupling.

The following figures illustrate these findings:

Figure 1: System-size scaling of structure factors, runtime, CG iterations per solve, and HMC acceptance rates for Holstein and oSSH models at $\beta$6.

Figure 2: Inverse-temperature scaling of structure factors, runtime, CG iterations, and HMC acceptance for $\beta$7 simulations.

Figure 3: Single-particle electron spectral function $\beta$8 for $\beta$9 lattices at $N$0, highlighting CDW gap formation (Holstein) and BOW-driven band renormalization (oSSH).

Statistical precision is achieved via long simulation runs and extensive averaging, while analytic continuation is performed with the DEAC algorithm, enabling precise extraction of low-energy spectral features relevant to experimental observables.

Implications, Comparisons, and Future Directions

Algorithmically, SmoQyElPhQMC.jl extends the current state of large-scale e-ph QMC by delivering the following practical advantages over legacy DQMC or pure HMC workflows:

• Orders-of-magnitude larger accessible lattice sizes,
• Adiabatic ($N$1) and strong-coupling regimes at finite doping,
• Efficient measurement infrastructure for momentum-resolved and real-space observables.

Theoretically, this package allows controlled studies of emergent phenomena in correlated electron-lattice systems, such as:

• Nature and criticality of CDW/BOW transitions in 2D,
• Nontrivial competition/coexistence between various symmetry-breaking ground states,
• Effects of dispersive phonons, nonlinearity, and disorder,
• Band structure renormalization and polaron formation at the thermodynamic limit.

Broader impacts are anticipated in:

• Materials theory, enabling direct model-based comparison to ARPES, RIXS, and transport experiments in quantum materials,
• Algorithmic development, as the preconditioning and pseudofermion integration strategies can be generalized to other finite-temperature lattice field theories (including gauge fields and multi-flavor Hubbard models) [Feng2025Scalable],
• Interoperability within the Julia quantum simulation ecosystem and extensibility to hybrid workflows (e.g., combining with machine-learning-accelerated Markov chain proposals [Albergo2021Flowbased]) and improved analytic continuation methods [Neuhaus2024SmoQyDEACjl].

Critical avenues for further research include:

• Systematic extension to correlated models with electron-electron interactions,
• Generalization of preconditioning, with potential use of Lanczos-based methods for further acceleration [Chen2024Lanczos],
• Automated tuning of HMC/EFA integrator parameters and deeper study of ergodicity solutions (radial update, sector crossing techniques) [Ostmeyer2025Exponential, Temmen2025Fully],
• Integration with high-level domain-specific libraries for automated model construction and experiment-theory comparison.

Conclusion

SmoQyElPhQMC.jl represents state-of-the-art methodology for QMC simulations of large-scale, uncorrelated e-ph Hamiltonians. Its innovative combination of HMC, exact Fourier acceleration, KPM preconditioning, and stochastic measurement infrastructure delivers near-linear scaling, high efficiency, model flexibility, and accurate physical observables—paving the way for systematic exploration of electron-lattice coupling in realistic materials contexts.

---

• "Fast and scalable quantum Monte Carlo simulations of electron-phonon models" [Cohen-Stead2022Fast]
• "SmoQyDQMC.jl: A flexible implementation of determinant quantum Monte Carlo for Hubbard and electron-phonon interactions" [Cohen-Stead2024SmoQyDQMCjl]
• "A hybrid Monte Carlo study of bond-stretching electron–phonon interactions and charge order in BaBiO$N$2" [Cohen-Stead2023Hybrid]
• "Charge order in the kagome lattice Holstein model: a hybrid Monte Carlo study" [Bradley2023Charge]
• "The Lanczos algorithm for matrix functions: a handbook for scientists" [Chen2024Lanczos]
• "Flow-based sampling for fermionic lattice field theories" [Albergo2021Flowbased]
• "SmoQyDEAC.jl: A differential evolution package for the analytic continuation of imaginary time correlation functions" [Neuhaus2024SmoQyDEACjl]

Additional relevant references provided in the main paper.