  We study the variational solution of generic interacting fermionic lattice
systems using fermionic Gaussian states and show that the process of
"gaussification", leading to a nonlinear closed equation of motion for the
covariance matrix, is locally optimal in time by relating it to the
time-dependent variational principle. By linearising our nonlinear equation of
motion around the ground-state fixed point we describe a method to study
low-lying excited states leading to a variational method to determine the
dispersion relations of generic interacting fermionic lattice systems. This
procedure is applied to study the attractive and repulsive Hubbard model on a
two-dimensional lattice.
