  The Lie-Trotter formula, together with its higher-order generalizations,
provides a direct approach to decomposing the exponential of a sum of
operators. Despite significant effort, the error scaling of such product
formulas remains poorly understood. We develop a theory of Trotter error that
overcomes the limitations of prior approaches based on truncating the
Baker-Campbell-Hausdorff expansion. Our analysis directly exploits the
commutativity of operator summands, producing tighter error bounds for both
real- and imaginary-time evolutions. Whereas previous work achieves similar
goals for systems with geometric locality or Lie-algebraic structure, our
approach holds in general. We give a host of improved algorithms for digital
quantum simulation and quantum Monte Carlo methods, including simulations of
second-quantized plane-wave electronic structure, $k$-local Hamiltonians,
rapidly decaying power-law interactions, clustered Hamiltonians, the transverse
field Ising model, and quantum ferromagnets, nearly matching or even
outperforming the best previous results. We obtain further speedups using the
fact that product formulas can preserve the locality of the simulated system.
Specifically, we show that local observables can be simulated with complexity
independent of the system size for power-law interacting systems, which implies
a Lieb-Robinson bound as a byproduct. Our analysis reproduces known tight
bounds for first- and second-order formulas. Our higher-order bound
overestimates the complexity of simulating a one-dimensional Heisenberg model
with an even-odd ordering of terms by only a factor of $5$, and is close to
tight for power-law interactions and other orderings of terms. This suggests
that our theory can accurately characterize Trotter error in terms of both
asymptotic scaling and constant prefactor.
