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.
This post has been migrated from (A51.SE)