BKT universality for height and planar spin fields

We prove Gaussian scaling limits for the two-dimensional discrete Gaussian height model throughout its rough phase, including the physical roughening threshold, for every finite square-symmetric interaction set containing the nearest neighbors. After the natural lattice normalization, the critical effective temperature has the universal value $8\pi$. For the ordinary nearest-neighbor Villain and XY models at sufficiently low fixed temperatures, we also prove that their Green-function-normalized spin fields converge to the imaginary exponential of a Dirichlet Gaussian free field.