Wilson loop operator in electrodynamics

1. Instead of showing Wilson loop is bounded, maybe showing its average is finite is more manipulable; Srednicki has a chapter discussing how to translate the average to a double-contour integral, which was first developed by Polyakov in 1980s.
2. The double-contour integral in the average naively is divergent, so combining the negative sign standing on the exponential function, the Mathematical average is zero. This is an intuitive estimation to your question.
3. Polyakov provided a procedure to regularize it. The subtraction term is known as perimeter law, and the regularized average is closely related to shape of the contour. But the value is bounded (from above) by double-contour integral around the circle.
4. Since $U(1)$ Wilson Loop reflects the physics of Schwinger effect, such boundness simply implies that the dynamic electric field provides more efficient production rate than static electric field. The background is also related to the technics, known as worldline instantons.

I think the approach you suggest in the question is probably doomed. Classically, the Wilson loop is the imaginary exponential $W_gamma(A) = e^{i int _gamma A}$. The map $xmapsto e^{ix}$ is bounded on the reals and holomorphic, so if you had a hermitian operator $A$, you'd be all set. But those operators don't exist in the quantum theory. (As an aside, if they did, there's no way they'd be bounded. They're going to get contributions from the large-field regions.)

What one can do rigorously is define the Maxwell theory via the generating function of expectation values of curvatures:

$$Z(h) = exp(-frac{1}{2} langle Delta^{-1} delta h, delta hrangle)$$

Here $Delta$ is the Laplacian on $1$-forms, $delta$ is the Hodge adjoint to the exterior derivative, and $h$ is a smooth, rapidly-decaying 2-form, the test function paired with the curvature $F_A$.

I'm going off memory here, so I've probably dropped some factors. But this generating function is a rigorous version of the generating functional $int e^{i langle A, dh rangle} e^{-frac{1}{4}||F_A||^2}dA$. You should think that $dh = J$ is the 3-form conserved current that naturally pairs with the connection $A$. Via the Bochner-Minlos theorem, it defines a Gaussian measure $dmu([A])$ on the linear dual of the space of smooth conserved currents $J$. The moments of this measure are Schwinger functions of the smeared curvature operators, e.g., $$int dx f(x)frac{d Z}{d h(x)} = E[F_A(f)],$$ for any reasonable test function $f in Omega^2$.

Osterwalder-Schrader applies, allowing one to interpet these Schwinger functions as operators on a Hilbert space. Then one can define the Wilson loop operator via area integral $W_gamma = exp(iint_R F_A)$ for some region $R$ bounded by your loop $gamma$. Then the reasoning from the first paragraph applies: This is a bounded function of a real (but probably unbounded) observable.