Optimized higher-order Lie-Trotter-Suzuki decompositions for two and more terms

Lie-Trotter-Suzuki decompositions of operator exponentials have a lot of applications in physics. For example, they are employed to sample equilibrium states in quantum Monte Carlo and to simulate the dynamics of quantum systems on quantum computers or on classical computers using tensor network state techniques. They also provide symplectic integrators for classical physics.
Good higher-order decompositions for exponentials of n=2 non-commuting operators are well-known. These cover one-dimensional quantum systems with nearest-neighbor interactions. We present some optimized decompositions for n=2 and new higher-order decompositions for n=3,4 which are needed for one-dimensional systems with longer-ranged interactions or quantum systems in higher dimensions. The ordering of operators in the decompositions turns out to have substantial influence on the attainable approximation order and magnitudes of the leading error terms.