Kevin Schäfers is a PhD student at the University of Wuppertal, in the Institute of Modelling,
Analysis, and Computational Mathematics (IMACM).
This post highlights the most recent results of the research project „Monte Carlo sampling of gluons by molecular dynamics“ (PI Michael Günther) within the research unit FOR5269 „Future methods for studying confined gluons in QCD“. The project primarily focuses on the design of efficient geometric numerical integration schemes for the Hamiltonian systems arising in the molecular dynamics step of the hybrid Monte Carlo (HMC) algorithm [7]. Within the HMC algorithm, the molecular dynamics step requires solving Hamiltonian equations of motion for a separable Hamiltonian H(p,q) = T(p) + V(q) with quadratic kinetic energy T(p).
The HMC algorithm
The HMC algorithm [7] is a frequent choice for gauge field generation in lattice QCD. The HMC algorithm consists of a molecular dynamics step and a Metropolis acceptance step.
The molecular dynamics step requires the solution of Hamiltonian equations of motion for a trajectory of length τ. In the standard HMC algorithm, the Hamiltonian is separable, H(p,q) = T(p) + V(q), and has a quadratic kinetic energy T(p). To satisfy the crucial detailed balance condition, ensuring that the Markov process converges to the correct equilibrium distribution, the numerical time integration schemes must be volume-preserving and time-reversible. Additionally, as the phase space of the Hamiltonian system is the cotangent bundle over a base space that is the special unitary group SU(3), numerical approximations must remain on this manifold.
Afterwards, in the Metropolis acceptance step, the new configuration is either accepted or rejected by investigating the difference in the Hamiltonian. Thus, another desirable property of the numerical integration scheme is good energy conservation, as this will result in a higher acceptance probability of the new configuration.
Splitting methods
For the numerical simulation of separable Hamiltonian systems, splitting methods provide a promising framework for geometric numerical integration [4,5]. Given the separability of the Hamiltonian, these methods can split the system into its kinetic and its potential part, each of which can be integrated exactly. Composing the flows of the respective subsystems yields a numerical approximation to the overall system that is symplectic (and thus volume-preserving) and, with a self-adjoint composition, time-reversible. Moreover, these integrators are explicit. These properties make splitting methods well-suited for the HMC’s molecular dynamics step.
For separable Hamiltonian systems with quadratic kinetic energy, force-gradient integrators [6] can provide even greater efficiency. Force-gradient integrators are a type of commutator-based splitting methods that incorporate an evaluation of a certain commutator into the computational process. As force-gradient integrators can be regarded as conventional splitting methods, applied to a Hamiltonian with modified potential, they are also frequently known as modified potential splitting methods. Moreover, this interpretation shows that force-gradient integrators do not compromise the geometric properties of splitting methods.
However, computing the required Hessian of the potential can be computationally costly or unavailable in software implementations. In such cases, Hessian-free force-gradient integrators [1] are a promising alternative as they approximate the force-gradient updates by using a second force evaluation, eliminating the occurrence of the Hessian and thus simplifying integration of these integrators into existing software packages. Due to the approximation, additional error terms occur and the integrators are no longer symplectic but still volume-preserving and hence still satisfy the crucial detailed balance condition.
Challenges in developing efficient integrators for lattice QCD simulations
The main goal in lattice QCD simulations is achieving a prescribed acceptance probability in the Metropolis step of the HMC algorithm at minimal computational cost. This requires tackling several challenges:
Higher convergence order. The total costs for ensemble generation scale with the lattice volume V = T x L3 as V4+1/2p with denoting the convergence order of the numerical integration scheme. Consequently, for increasing lattice sizes, higher-order integrators become more efficient. For instance, on a comparatively small lattice of size 48 x 243, fourth-order integrators turn out to be the most efficient choice [1-3].
Numerical stability. To achieve the target acceptance rate that is typically around 90%, we only require a coarse numerical approximation. However, the proposed splitting methods and (Hessian-free) force-gradient integrators are only conditionally stable. In the region of interest, the integrators are frequently affected by instabilities [3].
Multiple time stepping. The potential energy comprises a gauge and a fermion part, differing in magnitude and computational cost. As the fermion part causes the main part of the computational costs in the molecular dynamics step, the use of multiple time stepping techniques [8] employs different step sizes and saves costly evaluations of the fermion part while employing a sufficiently small step size to the gauge part. As a consequence, simulations in lattice QCD demand numerical time integration schemes at the intersection of geometric numerical integration and multirate integration.
Further splitting of the fermion part. Techniques like Hasenbusch mass preconditioning [9] split the fermion part further, each partition evaluated at different step sizes to avoid instabilities and improve efficiency.
All in all, we aim to design higher-order, Hessian-free force-gradient integrators that provide the optimal trade-off between accuracy and numerical stability, and incorporate multiple time
stepping techniques to avoid unnecessary evaluations of the costly fermion part.
The trade-off between accuracy and numerical stability
A common way to detect promising integrator variants within a certain order of convergence is to minimize a norm of the leading error coefficients, resulting in so-called minimum-error variants [6]. In [1], we used this approach to derive promising integrator variants within the framework of Hessian-free force-gradient integrators. When performing numerical tests, the asymptotic behavior matches the theoretical findings well. However, the approach assumes that all elementary differentials in the error expansion contribute equally to the overall discretization error. If the error terms are vastly different in magnitude, the chosen efficiency measure may not align with the experimental findings. In a next step, we will thus work on an integrator-based estimation of the error terms and then consider a problem-specific weighted norm of the leading error coefficients instead.
Next to the accuracy, the numerical stability of the integrators is of importance. By extending a linear stability analysis for splitting methods [10] to the framework of (Hessian-free) force-gradient integrators, we introduced a stability-related efficiency measure [3].
We found that the accuracy-related efficiency measure is way more sensitive to alterations in the integrator coefficients compared to the stability-related measure. We thus propose to keep maximizing the accuracy of a certain integrator variant by minimizing the norm of the leading error coefficients. Promising integrator variants are then obtained by selecting those minimum-error variants that also have a reasonably large stability threshold. Starting from 34 different variants of fourth-order Hessian-free force-gradient integrators listed in [1], we were able to find six non-dominated integrator variants, i.e., variants that are not outperformed by another integrator variant in both efficiency measures at the same time.
Outlook: extension to multiple time stepping
In lattice QCD, a common approach for multiple time stepping is nested integration. Nested integrators apply splitting methods in an iterative way and naturally incorporate an approach for multiple time stepping.
In a current work in progress, we formulate the idea of nested integration into a more general framework of hierarchical splitting methods that includes enhancements of the nested integration approach which are of particular interest in the derivation of higher-order multiple time stepping schemes. The framework also allows for making use of the force-gradient approach, as well as the Hessian-free variants. We moreover aim at generalizing the efficiency measures to the framework of hierarchical splitting methods.
References
[1] K. Schäfers, J. Finkenrath, M. Günther, F. Knechtli, Hessian-free force-gradient integrators, Computer Physics Communications 309 (2025), 109478. https://doi.org/10.1016/j.cpc.2024.109478
[2] K. Schäfers, J. Finkenrath, M. Günther, F. Knechtli, Hessian-free force-gradient integrators and their application to lattice QCD simulations, PoS LATTICE2024 (2025), 025. https://doi.org/10.22323/1.466.0025
[3] K. Schäfers, J. Finkenrath, M. Günther, F. Knechtli, Numerical stability of force-gradient integrators and their Hessian-free variants in lattice QCD simulations, arXiv preprint (2025). https://doi.org/10.48550/arXiv.2506.08813
[4] R.I. McLachlan, D.R.W. Quispel, Splitting Methods, Acta Numerica 11 (2002), pp. 341-434. https://doi.org/10.1017/S0962492902000053
[5] S. Blanes, F. Casas, A. Murua, Splitting methods for differential equations, Acta Numerica 33 (2024), pp. 1-161. https://doi.org/10.1017/S0962492923000077
[6] I.P. Omelyan, I.M. Mryglod, R. Folk, Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations, Computer Physics Communications 151.3 (2003), pp. 272-314. https://doi.org/10.1016/S0010-4655(02)00754-3
[7] S. Duane, A.D. Kennedy, B.J. Pendleton, D. Roweth, Hybrid Monte Carlo, Physics Letters B 195.2 (1987), pp. 216-222. https://doi.org/10.1016/0370-2693(87)91197-X
[8] C. Urbach, K. Jansen, A. Shindler, U. Wenger, HMC algorithm with multiple time scale integration and mass preconditioning, Computer Physics Communications 174.2 (2006), pp. 87-98. https://doi.org/10.1016/j.cpc.2005.08.006
[9] M. Hasenbusch, Speeding up the hybrid Monte Carlo algorithm for dynamical fermions, Physics Letters B 519.1-2 (2001), pp. 177-182. https://doi.org/10.1016/S0370-2693(01)01102-9
[10] S. Blanes, F. Casas, A. Murua, On the linear stability of splitting methods, Foundations of Computational Mathematics 8 (2008), pp. 357-393. https://doi.org/10.1007/s10208-007-9007-8