  We systematically extend Bogoliubov theory beyond the mean field
approximation of the Bose-Hubbard model in the superfluid phase. Our approach
is based on the time dependent variational principle applied to the family of
all Gaussian states (i.e. Gaussian TDVP). First, we find the best ground state
approximation within our variational class using imaginary time evolution in
1d, 2d and 3d. We benchmark our results by comparing to Bogoliubov theory and
DMRG in 1d. Second, we compute the approximate 1- and 2-particle excitation
spectrum as eigenvalues of the linearized projected equations of motion
(linearized TDVP). We find the gapless Goldstone mode, a continuum of
2-particle excitations and a doublon mode. We discuss the relation of the gap
between Goldstone mode and 2-particle continuum to the excitation energy of the
Higgs mode. Third, we compute linear response functions for perturbations
describing density variation and lattice modulation and discuss their relations
to experiment. Our methods can be applied to any perturbations that are linear
or quadratic in creation/annihilation operators. Finally, we provide a
comprehensive overview how our results are related to well-known methods, such
as traditional Bogoliubov theory and random phase approximation.
