  We derive the self-energy functional theory for bosonic lattice systems with
broken $U(1)$ symmetry by parametrizing the bosonic Baym-Kadanoff effective
action in terms of one- and two-point self-energies. The formalism goes beyond
other approximate methods such as the pseudoparticle variational cluster
approximation, the cluster composite boson mapping, and the Bogoliubov+U
theory. It simplifies to bosonic dynamical-mean field theory when constraining
to local fields, whereas when neglecting kinetic contributions of non-condensed
bosons it reduces to the static mean-field approximation. To benchmark the
theory we study the Bose-Hubbard model on the two- and three-dimensional cubic
lattice, comparing with exact results from path integral quantum Monte Carlo.
We also study the frustrated square lattice with next-nearest neighbor hopping,
which is beyond the reach of Monte Carlo simulations. A reference system
comprising a single bosonic state, corresponding to three variational
parameters, is sufficient to quantitatively describe phase-boundaries, and
thermodynamical observables, while qualitatively capturing the spectral
functions, as well as the enhancement of kinetic fluctuations in the frustrated
case. On the basis of these findings we propose self-energy functional theory
as the omnibus framework for treating bosonic lattice models, in particular, in
cases where path integral quantum Monte Carlo methods suffer from severe sign
problems (e.g. in the presence of non-trivial gauge fields or frustration).
Self-energy functional theory enables the construction of diagrammatically
sound approximations that are quantitatively precise and controlled in the
number of optimization parameters, but nevertheless remain computable by modest
means.
