Temperature-dependence of quark potential in Abelian lattice gauge theory

I am working with Kapusta's "Finite-Temperature Field Theory" textbook, and am working through the first part of chapter 10. When building the correlator of the two quarks a distance $$R$$ apart in the lattice gauge theory, he builds the partition function and finds

\begin{align} Z(\beta,\,R)&=\sum_\psi \langle \psi|e^{-\beta H}\int_{-\pi}^\pi \frac{d\alpha(x)}{2\pi} e^{i\alpha(0)}e^{-i\alpha(R)}|\psi\rangle\notag\\ &=\prod_x \int_{-\pi}^\pi \frac{d\alpha(x)}{2\pi} \prod_{{\textrm{links}}} \left( \sum_E \exp\left\{ -\frac{\beta^2 g}{2a}E^2+i (\alpha(x)-\alpha(x+n))E+i(\alpha(0)-\alpha(R))\right\} \right)\notag\\ &=Z(\beta)\langle e^{i\alpha(0)}e^{-i\alpha(R)}\rangle \notag\\ \end{align}

The free energy of the configuration is directly related to the thermal average contained in the above. In the low-temperature limit, he states that the flux can take on values of $$E=0$$ or $$\pm 1$$, and thus

\begin{align} &\langle e^{i\alpha(0)}e^{-i\alpha(R)}\rangle \notag\\ &=\prod_x \int_{-\pi}^\pi \frac{d\alpha(x)}{2\pi} \prod_{{\textrm{links}}} \left( \sum_E \exp\left\{ -\frac{\beta^2 g}{2a}E^2+i (\alpha(x)-\alpha(x+n))E+i(\alpha(0)-\alpha(R))\right\} \right)\notag\\ &=\prod_x \int_{-\pi}^\pi \frac{d\alpha(x)}{2\pi} \prod_{\textrm{links}}\left( e^{i(\alpha(0)-\alpha(R))}+e^{-\frac{\beta^2 g}{2a}+i\left\{\alpha(x)-\alpha(x+n)+\alpha(0)-\alpha(R)\right\}}+e^{-\frac{\beta^2 g}{2a}+i\left\{-\alpha(x)+\alpha(x+n)+\alpha(0)-\alpha(R)\right\}} \right) \end{align}

The first term is zero, while the second and third terms lead to a value of

$$\langle e^{i\alpha(0)}e^{-i\alpha(R)}\rangle =2e^{-\beta g^2 R/2a^2}$$

This I understand. However, in the high-temperature limit, we just have

$$\langle e^{i\alpha(0)}e^{-i\alpha(R)}\rangle =e^{\beta g^2/R}$$

I looked in this paper by Polyakov, which states that

$$\alpha(0)-\alpha(R)\sim -\frac{\beta g^2}{R}$$

and I'm not exactly sure why.

Ultimately, my question is how, exactly, do we get the coulombic potential at high temperature? In particular, why does Polyakov make the argument given above? Why can't we take $$E=0,\,\pm 1$$ as we did in the low-temperature case? Any explanation or references at the level of Kapusta would be appreciated.

EDIT: Also, any reference beyond those mentioned in the text that even talk about calculating the partition function for a lattice gauge theory in the high and low temperature regimes would be appreciated (online references are preferred).

A good reference which does the same calculation (at low temperature) in full QCD in a different way (which is more rigorous) is found in H.J. Rothe "Lattice Gauge Theory an introduction" in chapter 20 (and 18, 19 for some context) and also provides a calculation at high temperatures. For details over the calculations I suggest the references provided in this book. A couple of good references which may be useful are Critical behavior at finite-temperature confinement transitions by B.Svetitsky and L.G.Yaffe where the Coulomb gas formalism is developed and Phase transitions in Abelian lattice gauge theories by T.Banks, R.Myerson and J.Kogut.

The reason it is difficult to find a Polyakov-like derivation is probably that at high temperatures QCD is deconfining and perturbative techniques are, mostly, a valid way to estimate the $$q\bar{q}$$ potential which is what most people do.

Regarding the question at hand, I think that the reason you can consider only $$E=0,\pm1$$ is that in the Low temperature limit the eigenvalues of $$E$$ which dominate the exponential are the smallest ones. Since you have that $$\frac{\beta g^2}{2 a}>>1$$ $$\exp\left(-\frac{\beta g^2}{2 a} E^2+i \Delta \alpha\ E \right)$$ the exponential is highly suppressed and only the smallest eigenvalues contribute. Where $$\Delta \alpha \equiv \alpha(x)-\alpha(x+n)$$ for compactness.

On the other hand in the opposite limit $$\frac{\beta g^2}{2 a}<<1$$ such an approximation cannot hold since all eigenvalues will contribute to the exponential as long as $$\frac{\beta g^2}{2 a} E^2 \leq 1$$ In this case the path integral in $$E$$ will be dominated by the saddle point: $$E = \frac{i \Delta \alpha}{(\frac{\beta g^2}{2 a})}$$

Now we can try to recover the Polyakov result. It is possible to substitute the saddle point leading contribution and obtain

$$Z(R,\beta) = \int \mathcal{D} \alpha\ \prod_{links}\sqrt{\frac{2 \pi a}{\beta g^2}}\ e^{-\frac{a}{2 \beta g^2}\Delta \alpha^2}\times e^{i \alpha(0)-i\alpha(R)}$$

where $$\mathcal{D} \alpha\equiv \prod_x \int_{-\pi}^{\pi}\frac{d \alpha(x)}{2\pi}$$. Now we can approximate to the leading lattice contribution (i.e. ignoring lattice artifacts $$O(a)$$) and see that in this case $$\alpha$$ has a long range correlation: $$\prod_{links}\sqrt{\frac{2 \pi a}{\beta g^2}}\ e^{-\frac{a}{2 \beta g^2}\Delta \alpha^2} \sim\exp\left(-\frac{1}{\beta g^2}\int d^3x\ (\nabla \alpha)^2\right)$$ where the integral connects the links between the charges. The gradient pops out since $$\frac{1}{2a}\Delta \alpha\sim \nabla \alpha + O(a)$$ Therefore you see that $$Z(R,\beta) \sim\int \mathcal{D} \alpha\ \exp\left(-\frac{1}{ \beta g^2}\int d^3x\ (\nabla \alpha)^2 \right)\ e^{i \alpha(0)-i\alpha(R)}$$ where now the integration generates the propagator connecting $$0$$ to $$R$$. To see this, just imagine that $$\alpha(0) = \int \alpha(x)J(x)$$ where $$J(x) = \delta(x)$$ and the same for $$\alpha(R)$$ then you can use standard QFT results about gaussian integrals with sources to obtain: $$Z(R,\beta) \sim \exp\left(-\beta g^2 G(0,R)\right)$$ where $$\nabla^2 G(0,R)=\delta(R)$$ and therefore in coordinate space $$Z(R,\beta) \sim \exp\left(-\frac{\beta g^2}{R}\right)$$ up to constant terms.