Consider the XY model on the square lattice. A field configuration $\theta$ is specified by an element of the Abelian group $\mathbb{R}/2\pi \simeq U(1)$ at each vertex of the lattice. The gradient of this field is $$(\mathrm{d}\theta)_{ij} := \theta_j - \theta_i,$$ where $\mathrm{d}$ is the lattice exterior derivative and $i,j$ are two vertices bounding a single edge. Note that the difference here is taken in $\mathbb{R}/2\pi$. Then for each plaquette $p$ with sequential corners $i,j,k,l$, we have $$(\mathrm{d}^2\theta)_{p} := (\theta_j-\theta_i) + (\theta_k-\theta_j) + (\theta_l - \theta_k) + (\theta_i - \theta_l) \in \mathbb{Z} \,\,(\text{mod } 2\pi)$$ Actually, it can only be $0,\pm 1$, if for each difference we take $\theta_j - \theta_i \in [-\pi,\pi)$, and the sum is taken in $\mathbb{R}$: the result is zero if we take all sums mod $2\pi$.
My primary question is, how exactly do I formalize the fact that this quantity is an integer mod $2\pi$ and not zero? My only thought is that $$(\mathrm{d}^2\theta)(p) = (\mathrm{d}\theta)(\partial p)$$ where $\partial$ denotes boundary and somehow this defines a map $S^1 \to U(1)$ and therefore a homotopy class, but I'm struggling to make this explicit. Non-trivial solutions of $\mathrm{d}^2 = 0$ in the continuum are only possible for singular field configurations, but these vortices are topological on the lattice in some well defined sense that I am not understanding.
As a follow up question, if I replace $\mathbb{R}/2\pi$ by $\mathbb{Z}/n\mathbb{Z}$ (a clock model) will I also get $\mathbb{Z}_n$ vortices? In that case there is no notion of a homotopically non-trivial map $S^1 \to \mathbb{Z}_n$, but perhaps there is a discrete generalization of homotopy that applies?
