  We present a new variational method for investigating the ground state and
out of equilibrium dynamics of quantum many-body bosonic and fermionic systems.
Our approach is based on constructing variational wavefunctions which extend
Gaussian states by including generalized canonical transformations between the
fields. The key advantage of such states compared to simple Gaussian states is
presence of non-factorizable correlations and the possibility of describing
states with strong entanglement between particles. In contrast to the commonly
used canonical transformations, such as the polaron or Lang-Firsov
transformations, we allow parameters of the transformations to be time
dependent, which extends their regions of applicability. We derive equations of
motion for the parameters characterizing the states both in real and imaginary
time using the differential structure of the variational manifold. The ground
state can be found by following the imaginary time evolution until it converges
to a steady state. Collective excitations in the system can be obtained by
linearizing the real-time equations of motion in the vicinity of the imaginary
time steady-state solution. Our formalism allows us not only to determine the
energy spectrum of quasiparticles and their lifetime, but to obtain the
complete spectral functions and to explore far out of equilibrium dynamics such
as coherent evolution following a quantum quench. We illustrate and benchmark
this framework with several examples: a single polaron in the Holstein and
Su-Schrieer-Heeger models, non-equilibrium dynamics in the spin-boson and Kondo
models, the superconducting to charge density wave phase transitions in the
Holstein model.
