  Variational minimization of tensor network states enables the exploration of
low energy states of lattice gauge theories. However, the exact numerical
evaluation of high-dimensional tensor network states remains challenging in
general. In [E. Zohar, J. I. Cirac, Phys. Rev. D 97, 034510 (2018)] it was
shown how, by combining gauged Gaussian projected entangled pair states with a
variational Monte Carlo procedure, it is possible to efficiently compute
physical observables. In this paper we demonstrate how this approach can be
used to investigate numerically the ground state of a lattice gauge theory.
More concretely, we explicitly carry out the variational Monte Carlo procedure
based on such contraction methods for a pure gauge Kogut-Susskind Hamiltonian
with a $\mathbb{Z}_3$ gauge field in two spatial dimensions. This is a first
proof of principle to the method, which provides an inherent way to increase
the number of variational parameters and can be readily extended to systems
with physical fermions.
