Numerical stability of force-gradient integrators and their Hessian-free variants in lattice QCD simulations

TL;DR

This paper conducts a detailed stability analysis of force-gradient and Hessian-free integrators in lattice QCD, identifying variants that optimize accuracy and stability for efficient simulations, especially in Hamiltonian Monte Carlo methods.

Contribution

It provides a comprehensive stability analysis of force-gradient integrators and their Hessian-free variants, and demonstrates their effectiveness in lattice QCD simulations.

Findings

Hessian-free and conventional force-gradient integrators share similar stability domains.

Certain integrator variants offer a good balance between accuracy and stability.

Hessian-free force-gradient integrators enable more efficient lattice QCD simulations.

Abstract

A comprehensive linear stability analysis of force-gradient integrators and their Hessian-free variants is carried out by investigating the harmonic oscillator as a test equation. The analysis reveals that the linear stability of conventional force-gradient integrators and their Hessian-free counterparts coincides. By performing detailed linear stability investigations for the entire family of self-adjoint integrators with up to eleven exponentials per time step, we detect promising integrator variants that are providing a good trade-off between accuracy and numerical stability. Special attention is given to the application of these promising integrator variants within the Hamiltonian Monte Carlo algorithm, particularly in the context of interacting field theories. Simulations for the two-dimensional Schwinger model are conducted to demonstrate that there are no significant differences in the stability domain of a force-gradient integrator and its Hessian-free counterpart. Lattice QCD simulations with two heavy Wilson fermions emphasize that Hessian-free force-gradient integrators with a larger stability threshold allow for a more efficient computational process compared to conventional splitting methods. Furthermore, detailed investigations of the stability threshold are performed by investigating Nf = 2 twisted-mass fermions and nested integrators, highlighting the reliability of the linear stability threshold for lattice QCD simulations.