Relation between sources and Wilson loops in Tong's gauge theory In Euclidean quantum Yang-Mills with compact gauge group $G$, the VEV of a Wilson loop is:
$$\tag{$\star$} \langle W[C]\rangle \equiv \int_{\mathcal{A}/G}DA e^{-S_{YM}[A]}tr\mathcal{P}\exp\left(i\int_C A\right).$$
On p87 of David Tong's notes on gauge theory (in the chapter entitled Yang-Mills Theory) he writes:
Now consider the specific closed loop C shown in the figure. We again take this to sit in
the fundamental representation. It has the interpretation that we create a quark antiquark pair, separated by a distance r, at some time in the past. These then propagate
forward for time T, before they annihilate back to the vacuum.

Then Tong says something crucial which I don't understand:
In the presence of the sources, the ground state of the system has energy $V(r)$.
My confusion isn't about the specific form of $V(r)$ or anything like that; I just don't know what the words in the above sentence are meant to mean.
If we were studying a scalar field, say, I'd have no issue. In that case, "in the presence of sources" would mean "insert some $\exp(\int J(x)\phi(x))$ into the path integral". For our specific scenario we'd take something like $J(x)=\delta(\vec{x})+\delta(\vec{x}-\vec{r})$, representing static sources at $\vec{0}$ and $\vec{r}$. Then "the ground state of the system has energy $V(r)$" means that if we include the sources in the Lagrangian, then the resulting Hamltonian $H_J$ has ground state energy $V(r)$. We can compute $V(r)$ perturbatively, with the leading order given in terms of the tree propagator.
In Yang-Mills, such interpretations are unavailable. Gauge invariant operator insertions are given by Wilson loops, not source terms of the form $\int J(x)\phi(x)$. Moreover, the Wilson loop only exists "as a whole"; it doesn't make sense to take a slice of time and ask what the Hamiltonian is in the presence of that Wilson loop. So I don't know what "ground state energy in the presence of sources" means here, let alone how to calculate it.
To summarise my question, what does Tong mean by the statement: In the presence of the sources, the ground state of the system has energy $V(r)$?
Ideally, I'd like an answer which works for strong coupling, since we'd like to conclude that an area law for the Wilson loop VEV implies a confining force $V(r)=\sigma r$: a fact which is supposed to be true in the strongly-coupled regime.
 A: Wilson lines indeed act as sources for charged particles. For the purposes of illustration, let me go to the abelian case momentarily, so that $\mathrm{tr}\,\mathrm{P}\!\exp(\cdot)$ is simply $\exp(\cdot)$, and I'll return to Yang-Mills in the end.
$\newcommand{\d}{\mathrm{d}}$
Abelian case
So consider an abelian gauge theory with a gauge field, $a$, and now you want to add a source. The obvious generalisation of "inserting some $\exp(\mathrm{i}\;J\cdot\phi)$ in the path integral" is inserting some
$$\exp\!\left(\mathrm{i}\;\int_\mathcal{M}\d^d x\ J^\mu(x)\ a_\mu(x)\right)\tag{1}$$
in the path integral. And now, let's come to differential forms. The gauge field, $a = a_\mu(x)\;\d x^\mu$, is a one-form and so is $J = J_\mu(x)\; \d x^\mu$. The integral in the expression (1) can be rewritten as
$$ \int_\mathcal{M}\d^d x\ J^\mu(x)\ a_\mu(x) = \int_\mathcal{M} a\wedge\star J =: \int_\mathcal{M} a\wedge K,$$
where $\star$ is the Hodge-star operator and I defined the $(d-1)$-form
$$K :=\star J = K_{\mu_1 \mu_2\cdots \mu_{d-1}} \d{x}^{\mu_1} \wedge \d{x}^{\mu_2} \wedge \cdots \wedge \d{x}^{\mu_{d-1}}.$$
Now we come to the fun part, Poincaré duality! Since your theory has gauge invariance, $a\sim a+\d \lambda $ an integration by parts shows that $\d\star J = 0$, or equivalently $\d K = 0$. That means that $K$ is a representative of a cohomology class of the $(d-1)$-th cohomology group of $\mathcal{M}$. Invoking Poincaré duality, i.e. $\mathrm{H}^n(\mathcal{M}) \cong \mathrm{H}_{d - n}(\mathcal{M})$ you can rewrite the above integral as
$$ \int_\mathcal{M} a\wedge K = \int_\widehat{K} a, \tag{2}$$
where $\widehat{K}$ is the Poincaré dual 1-cycle to the $(d-1)$-class $[K]$ (see appendix for a proof of (2)). That's it! If I rename $\widehat{K}$ as $C$, this argument just shows that
$$\exp\!\left(\mathrm{i} \int_C a\right) =: W[C],$$
is a source for charged particles.
In the opposite direction, starting with a Wilson loop $W[C]$ you can reverse engineer $J$ from (1), as
$$J_\mu(x) = (-1)^{(d-1)+\epsilon}\big[\star\!\widehat{C}\,\big]_\mu(x),$$
where $\epsilon=0$ if you are in Euclidean signature and $\epsilon=1$ in Lorentzian, and $\widehat{C}$ is (a representative of) the Poincaré dual $(d-1)$-class, of the 1-cycle $C$. This is roughly speaking a $(d-1)$-form delta function, restricting you to lie in $C:\quad$ $ \widehat{C} = \delta_{[d-1]}\big(x\in C\big), $ so
$$ J_\mu(x) = (-1)^{(d-1)+\epsilon} \left[\star\delta_{[d-1]}\big(x\in C\big)\right]_\mu.$$
Back to Yang-Mills
All of the above differential-geometric arguments hold the same in a non-abelian gauge theory, by replacing exponentials with traces of path-ordered exponentials, therefore the slogan is inserting a Wilson loop
$$W[C]:= \mathrm{tr}\,\mathrm{P}\!\exp\!\left(\mathrm{i}\int_C a\right)$$
in the path-integral, is equivalent to adding a source given by a current:
$ J_\mu(x) = (-1)^{(d-1)+\epsilon} \left[\star\delta_{[d-1]}\big(x\in C\big)\right]_\mu.$ If you wish to find the ground state energy in the presence of a Wilson loop, everything is now ready to be translated into Hamiltonian language, exactly like you did for the scalar case.
The path-ordering symbol somewhat obscures the source interpretation. One way to get around this, is described also in Tong's gauge theory lecture notes, around subsection 2.1.3. The trick is to trade the path-ordering symbol for auxiliary fields. In particular, you can write:
\begin{align} W_\mathrm{R}[C] &= \mathrm{tr}_\mathrm{R}\,\mathrm{P}\!\exp\!\left(\mathrm{i}\int_C a\right) \\
&= \int \mathrm{D}\lambda \mathrm{D}w^\dagger \mathrm{D}w \ \mathcal{O}_\mathrm{R}\left(w,w^\dagger\right)\ \exp\!\left(\mathrm{i}\int_C \left(w^\dagger \mathrm{i}\,\d w + \lambda (w^\dagger w-\kappa_\mathrm{R}) + w^\dagger w\; a\right) \right), \tag{3}
\end{align}
where $w$ is a complex vector of dimension $\dim\mathrm{R}$ ($\mathrm{R}$ is the representation on which the Wilson loop is computed), $\kappa_\mathrm{R}$ is number that depends on the representation, that controls the length of $w$, and $\mathcal{O}_\mathrm{R}$ is an insertion whose exact form depends on the representation. From here on, you can play the Poincaré game, make $C$ into a $J$ and you are left with a source term [the last term in (3)] and a bunch of auxiliary fields that you need to integrate out.

Appendix: Justifying Poincaré duality
$\newcommand{\cM}{\mathcal{M}}$In this appendix I will justify equation (2), namely, I will show that Poincaré duality can be used with $a$ being any $1$-form. Let's define the pairing
\begin{align}
\Omega^1(\cM)\times\Omega^{d-1}(\cM)&\to \mathbb{R} \\ 
\left(\omega,\eta\right)&\mapsto \int_\cM \omega\wedge\eta.
\end{align}
Fixing an element $a\in\Omega^1(\cM)$ and restricting the second factor to a cohomology class, induces a map
\begin{align}
\phi_a:\mathrm{H}^{d-1}(\cM) &\to \mathbb{R} \\
\eta &\mapsto \phi_a(\eta) = \int_\cM a\wedge \eta.
\end{align}
I can also construct a map
\begin{align}\psi_a :\mathrm{H}_{1}(\cM) &\to \mathbb{R} \\
C&\mapsto \psi_a(C):=\int_C a,
\end{align}
which is simply integrating $a$ on a 1-cycle. Poincaré duality asserts that $$\mathrm{H}_{1}(\cM) \cong \mathrm{H}^{d-1}(\cM).$$
I can rewrite that by saying that there is an isomorphism
$$\xi :\mathrm{H}_{1}(\cM) \overset{\sim}{\longrightarrow} \mathrm{H}^{d-1}(\cM).$$
Finally I can construct the following commutative square of maps:
$$\require{AMScd}
\begin{CD}
\mathrm{H}_{1}(\cM) @>{\xi}>{\sim}> \mathrm{H}^{d-1}(\cM) \\
@V{\psi_a}VV @VV{\phi_a}V \\
\mathbb{R}@>{\sim}>{\mathrm{id}}>\mathbb{R},
\end{CD}$$
which shows that $\phi_a\circ\xi = \mathrm{id}\circ\psi_a=\psi_a$, or in other words, for any $a\in\Omega^1(\cM)$
$$\int_C a=:\psi_a(C) = \phi_a(\xi(C)) := \int_\cM a\wedge \xi(C).$$
This is precisely equation (2) with $\xi(\widehat{K})=K$.
