Given a class of $q$-local Hamiltonians, is it possible to find a simple variational state whose energy is a finite fraction of the ground state energy in the thermodynamic limit? Whereas product states often provide an affirmative answer in the case
of bosonic (or qubit) models, we show that Gaussian states fail dramatically in the fermionic case, like for the Sachdev-Ye-Kitaev (SYK) models. This prompts us to propose a new class of wavefunctions for SYK models inspired by the variational coupled cluster algorithm. We introduce a static (0+0D) large-$N$ field theory to study the energy, two-point correlators, and entanglement properties of these states. Most importantly, we demonstrate a finite disorder-averaged approximation ratio of $r approx 0.62$ between the variational and ground state energy of SYK for $q=4$. Moreover, the variational states provide an exact description of spontaneous symmetry breaking in a related two-flavor SYK model.
