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!