  We introduce the harmonic virtual element method (harmonic VEM), a
modification of the virtual element method (VEM) for the approximation of the
2D Laplace equation using polygonal meshes. The main difference between the
harmonic VEM and the VEM is that in the former method only boundary degrees of
freedom are employed. Such degrees of freedom suffice for the construction of a
proper energy projector on (piecewise harmonic) polynomial spaces. The harmonic
VEM can also be regarded as an "$H^1$-conformisation" of the Trefftz
discontinuous Galerkin-finite element method (TDG-FEM). We address the
stabilization of the proposed method and develop an $hp$ version of harmonic
VEM for the Laplace equation on polygonal domains. As in Trefftz DG-FEM, the
asymptotic convergence rate of harmonic VEM is exponential and reaches order
$\mathcal O ( \exp(-b\sqrt[2]{N}))$, where $N$ is the number of degrees of
freedom. This result overperformes its counterparts in the framework of $hp$
FEM and $hp$ VEM, where the asymptotic rate of convergence is of order
$\mathcal O ( \exp(-b\sqrt[3]{N}) )$.
