Expansion of Wilson loop operator to obtain the field tensor in non-Abelian gauge theory

Question

I am reading Prof. Eduardo Fradkin's QFT lecture note (the link to the lecture note is gone). In his note, he considered the Wilson loop operator
\begin{equation}
\hat{W}_{\Gamma(x,x)}=\hat{P} \left( \exp{ig\oint_{\Gamma(x,x)}dz^{\mu}A_{\mu}(z)} \right)
\end{equation}
where $\Gamma(x,x)$ denotes a loop with both its initial and final points at $x$, $g$ is the coupling constant, and $A_{\mu}$ is the non-Abelian gauge field. We want to obtain the definition of the field tensor $F_{\mu \nu}$ by expanding $W_{\Gamma(x,x)}$ to the leading order when we have infinitesimal $\Gamma(x,x)$. This means the length of the loop $\Gamma(x,x)$ is infinitesimal and also the oriented surface whose boundary is $\Gamma(x,x)$ has its minimal area infinitesimal. We then proceed as follows
\begin{equation}
\hat{W}_{\Gamma(x,x)} \approx \hat{P} \left( 1 + ig\oint_{\Gamma(x,x)}dz^{\mu}A_{\mu}(z) + \frac{1}{2} \left[ ig\oint_{\Gamma(x,x)}dz^{\mu}A_{\mu}(z) \right]^{2} + \cdots \right).
\end{equation}
The first two terms are easy to do via Stoke's theorem:
\begin{align}
& \hat{P} \left( 1\right) = 1 \\
& \hat{P} \left( ig\oint_{\Gamma(x,x)}dz^{\mu}A_{\mu}(z) \right) = ig\oint_{\Gamma(x,x)}dz^{\mu}A_{\mu}(z) = \frac{ig}{2} \int_{\Sigma} dx^{\mu} \wedge dx^{\nu} \left( \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu} \right).
\end{align}
My question is that he then directly wrote down
\begin{align}
\hat{P} \left( \frac{1}{2}  \left[ \oint_{\Gamma(x,x)}dz^{\mu}A_{\mu}(z) \right]^{2} \right) = \frac{1}{2} \int_{\Sigma} dx^{\mu} \wedge dx^{\nu} \left( -\left[A_{\mu},A_{\nu} \right] \right).
\end{align}
**I am a bit confused about how can we turn the path-ordered product of two line integrals into something like a surface integral.** I currently do not have an idea about where to start to prove this identity. But it seems to be crucial if we want to prove that, for infinitesimal loop $\Gamma(x,x)$, we have
\begin{align}
\hat{W}_{\Gamma(x,x)} \approx 1 + \frac{ig}{2}\int_{\Sigma} dx^{\mu} \wedge dx^{\nu} F_{\mu\nu}
\end{align}
where
\begin{align}
F_{\mu\nu} = \partial_{\mu} A_{\nu} - \partial_{\nu} A_{\mu} -ig[A_{\mu},A_{\nu}].
\end{align}
Are there any suggestions or good reference that carry out this computation
\begin{align}
\hat{P} \left( \frac{1}{2}  \left[ \oint_{\Gamma(x,x)}dz^{\mu}A_{\mu}(z) \right]^{2} \right) = \frac{1}{2} \int_{\Sigma} dx^{\mu} \wedge dx^{\nu} \left( -\left[A_{\mu},A_{\nu} \right] \right)
\end{align}
in details? Thanks!