Derivation of the effective potential between a quark and an anti-quark Typically in particle physics books (not in QFT books!) I have often seen  this statement that the potential between a heavy quark and its anti-quark can be "empirically" represented as $V(r) = -\frac{\alpha_s}{r} + br$ where $\alpha_s \sim \frac{\hbar c}{2}$ and $b \sim \frac{0.18}{\hbar c} GeV^2$.


*

*Is there a way to get the above form or the expression from a QCD calculation?


I have seen some approximate evaluations of the Wilson loops whereby one shows that the $2$-point gauge field correlator decays as $\frac{1}{r}$ to quadratic order of the coupling but thats about it. And I am not sure if somehow the proof that the beta-function of $QCD$ is decreasing to first order somehow implies the above form for $V(r)$. 
I would like to know how to derive the above expression for $V(r)$. 
 A: Recent developments on lattice computations in QCD have shown that the beta function for a pure Yang-Mills theory (let me emphasize that is not true for the full QCD) goes to zero lowering momenta and so, the theory seems to reach a trivial infrared fixed point (see here, fig. 5). The matter about low-energy behavior of a Yang-Mills theory seems rather well settled through studies as the one I have just cited.
A trivial infrared fixed point means that the theory is well described by a free one at that limit and, if one knows the form of the propagator, a form of the potential can be computed through Wilson loop (see here). 
I give this computation here. For a free theory, the case of the infrared trivial fixed point, the generating functional is Gaussian (I will consider Landau gauge)
\begin{equation}
     Z_0[j]=\exp\left[\frac{i}{2}\int d^4x'd^4x''j^{a\mu}(x')D_{\mu\nu}^{ab}(x'-x'')j^{b\nu}(x'')\right]
\end{equation}
and so, the evaluation of the Wilson loop is straightforwardly obtained as
\begin{equation}
     W[{\cal C}]=\exp\left[-\frac{g^2}{2}C_2\oint_{\cal C}dx^\mu\oint_{\cal C}dy^\nu D_{\mu\nu}(x-y)\right]
\end{equation}
being $C_2=$ the quadratic Casimir operator that for SU(N) is $(N^2-1)/2N$. The fall-off to large distances of this propagator grants that ordinary arguments to evaluate the integrals on the path apply. Indeed, using Fourier transform one has
\begin{equation}
     W[{\cal C}]=\exp\left[-\frac{g^2}{2}C_2\int\frac{d^4p}{(2\pi)^4}\Delta(p)\left(\eta_{\mu\nu}-\frac{p_\mu p_\nu}{p^2}\right)\oint_{\cal C}dx^\mu\oint_{\cal C}dy^\nu e^{-ip(x-y)}\right].
\end{equation}
We need to evaluate
\begin{equation}
     I({\cal C})=\eta_{\mu\nu}\oint_{\cal C}dx^\mu\oint_{\cal C}dy^\nu e^{-ip(x-y)}
\end{equation}
provided the contributions coming from taking into account the term $\frac{p_\mu p_\nu}{p^2}$ run faster to zero at large distances. This must be so also in view of the gauge invariance of Wilson loop. With the choice of time component in the loop going to infinity while distance is kept finite, we can evaluate the above integral in the form
\begin{equation}
     I({\cal C})\approx 2\pi T\delta(p_0)e^{-ipx}
\end{equation}
and we are left with
\begin{equation}
     W[{\cal C}]\approx \exp\left[-T\frac{g^2}{2}C_2\int\frac{d^3p}{(2\pi)^3}\Delta({\bf p} ,0)e^{-i{\bf p}\cdot{\bf x}}\right]
\end{equation}
yielding
\begin{equation}
     W[{\cal C}]=\exp\left[-TV_{YM}(R)\right].
\end{equation}
So, if you are able to compute Yang-Mills propagator in the low-energy limit you will get an answer to your problem.
