  The Luttinger-Ward functional $\Phi[\mathbf{G}]$, which expresses the
thermodynamic grand potential in terms of the interacting single-particle
Green's function $\mathbf{G}$, is found to be ill-defined for fermionic models
with the Hubbard on-site interaction. In particular, we show that the
self-energy $\mathbf{\Sigma}[\mathbf{G}] \propto \delta\Phi[\mathbf{G}]/\delta
\mathbf{G}$ is not a single-valued functional of $\mathbf{G}$: in addition to
the physical solution for $\mathbf{\Sigma}[\mathbf{G}]$, there exists at least
one qualitatively distinct unphysical branch. This result is demonstrated for
several models: the Hubbard atom, the Anderson impurity model, and the full
two-dimensional Hubbard model. Despite this pathology, the skeleton Feynman
diagrammatic series for $\mathbf{\Sigma}$ in terms of $\mathbf{G}$ is found to
converge at least for moderately low temperatures. However, at strong
interactions, its convergence is to the unphysical branch. This reveals a new
scenario of breaking down of diagrammatic expansions. In contrast, the bare
series in terms of the non-interacting Green's function $\mathbf{G}_0$
converges to the correct physical branch of $\mathbf{\Sigma}$ in all cases
currently accessible by diagrammatic Monte Carlo. Besides their conceptual
importance, these observations have important implications for techniques based
on the explicit summation of diagrammatic series.
