  The time-evolving matrix product operators (TEMPO) method, which makes full
use of the Feynman-Vernon influence functional, is the state-of-the-art tensor
network method for bosonic impurity problems. However, for fermionic impurity
problems the Grassmann path integral prohibits application of this method. We
develop Grassmann time-evolving matrix product operators, a full fermionic
analog of TEMPO, that can directly manipulates Grassmann path integrals with
similar numerical cost as the bosonic counterpart. We further propose a zipup
algorithm to compute expectation values on the fly without explicitly building
a single large augmented density tensor, which boosts our efficiency on top of
the vanilla TEMPO. Our method has a favorable complexity scaling over existing
tensor network methods, and we demonstrate its performance on the
non-equilibrium dynamics of the single impurity Anderson models. Our method
solves the long standing problem of turning Grassmann path integrals into
efficient numerical algorithms, which could significantly change the
application landscape of tensor network based impurity solvers, and could also
be applied for broader problems in open quantum physics and condensed matter
physics.
