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

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.