Quantization of Flux in Polyakov's 3D Compact QED In his his book "Gauge Fields and Strings" Polyakov introduces the compact QED on a cubic lattice in 3D Euclidean space as:
$$ S\left[ \left\{ A_{\mathbf{r},\mathbf{\alpha}}\right\} \right]=\frac{1}{2g^2}\sum_{\mathbf{r},\mathbf{\alpha},\mathbf{\beta}}(1-\cos{F_{\mathbf{r},\mathbf{\alpha}\mathbf{\beta}}}) $$
Where $F$ is the net flux through the plaquette which is spanned by the lattice vectors $\mathbf{\alpha}$ and $\beta$ at point $\mathbf{r}$ and is given by:
$$ F_{\mathbf{r},\mathbf{\alpha}\mathbf{\beta}}=A_{r,\alpha}+A_{r+\alpha,\beta}-A_{r,\beta}-A_{r+\beta,\alpha}$$
Which intuitively is the curl of $A$ around the plaquette. The gauge transformation is defined as:
$$ A_{r,\alpha}\to A_{r,\alpha}-\phi_{r}+\phi_{r+\alpha} $$
Under which the action is invariant.
One obvious result is that the total flux through any closed Gaussian surface is zero. This is true because:
$$\sum_{p\in cube} F_p=0$$
As each gauge field on each link appears twice with different signs in the above sum. So it is impossible to have monopoles in this system except for Dirac monopoles which can be built by assuming that the flux through 5 faces of a cube has the same sign while one face has a net flux with negative sign such that the total flux remains zero.
But then, he (Polyakov) states that this flux (which only goes through one of the faces of a cube) is quantized. I do not know how to prove this. It seems that a singular gauge transformation is necessary (according to a paper by 't Hooft) and we need to couple the gauge field to another (probably matter) field, but I cannot find a way to implement that transformation in the lattice model and even one might ask why should we couple $A$ to another degree of freedom. This point is also mentioned here: https://physics.stackexchange.com/a/202806/90744 again without any proof.
The book uses another action which is claimed to be equivalent to the original action, which is given by:
$$ S=\frac{1}{4g^2}\sum_{r,\alpha,\beta}(F_{r,\alpha \beta}- 2\pi n_{r,\alpha \beta})^2 $$
Where $n$ is an integer valued field. This action in general is not equivalent to the original action. because here we are allowing deviations from non-periodicity of $A$ to contribute and therefore we can only use it in the small $g$ limit.
 A: Well, concerning the question, it should follow from the discrete version of the Stokes' theorem. Consider a cube, in case of a nonzero flux, piercing the cube, one cannot assign globally the gauge potential $A_\mu$, only locally, in a certain chart. Let us divide the cube into two charts, overlapping at least on the equator 
The northern and southern hemisphere. According to the Stokes' theorem flux through the pale red surface is equal to the circulation of $A_\mu$ around the equator:
$$
\int_{U_N} F d S= \sum_{i \in s} F_i S_i = \oint A_\mu dx^{\mu} = \sum_{i \in l} A_i l_i 
$$
Where $s$ - denotes all surfaces in the chart, and $l$ - the line segments on the equator, and $S_i$ - area of the surface, $l_i$ - length of the segment. In the integral over equator, one may choose in Stokes' theorem to integrate over $U_N$ and $U_S$, and the result, from the physical point of view, should not depend on the choice of surface.
The electromagnetic part of the action for the point particle is:
$$
S = \oint A_\mu d x^{\mu}
$$
The action for the point particle enters the path integral as $e^{i S}$
Therefore, in order for the $e^{i S}$ to be single-valued, the fluxes over northern and southern hemisphere have to satisfy following condition:
$$
\int _{U_N} F = - \int _{U_S} F + 2 \pi n \qquad n \in \mathbb{Z} \qquad
\Rightarrow \qquad \int _{U_N \cup U_S} F = 2 \pi n 
$$
This logic lacks rigor, but may provide some intuition. Another point, which one can note, that the monopoles are the classical solutions - minima of the action functional, and from the action, one can see, that:
$$
\cos F_{r, \alpha \beta} = 1 \Rightarrow F_{r, \alpha \beta} = 2 \pi n, n \in \mathbb{Z}
$$
So the sum over all faces, will be quantized.
The action, that you have written in the end of your post, is a Villain or Gaussian approximation of the original action, which assumes, that fluctuations of the gauge field are close to the minima $F_{r, \alpha \beta} = 2 \pi n$, and is obtained by the expansion of the cosine to the second order:
$$
1 - \cos F_{r, \alpha \beta} = \frac{1}{2} (F_{r, \alpha \beta} - 2 \pi n_{r, \alpha \beta})^2
$$
