Intersecting Wilson loops in 2D Yang-Mills I am currently trying to understand 2D Yang-Mills theory, and I cannot seem to find an explanation for calculation of the expectation value of intersecting Wilson loops.
In his On Quantum Gauge Theories in two dimensions, Witten carries out a curious calculation:
For three reps $ \alpha,\beta,\gamma $, we fix a basis of the tensor product space belonging to $\alpha \otimes \beta \otimes \gamma$ called $\epsilon_\mu(\alpha\beta\gamma)^{ijk}$ ($\mu$ indexes the $\mu$-th basis vector, the $i,j,k$ are the indices from the original reps) with the property that
$$ \int \mathrm{d}U {\alpha(U)^i}_i' {\beta(U)^j}_j' {\gamma(U)^k}_k' = \epsilon_\mu(\alpha\beta\gamma)^{ijk}\bar{\epsilon}_\mu(\alpha\beta\gamma)_{i'j'k'}$$
A minor question is why this is possible - I would be fine with accepting that I can always find some vectors that fulfill that relation, but why are they a basis?
The real part I do not understand comes now: By the above, each edge of a plaquette carries some $\epsilon_\mu$, and at a crossing of two lines, we have thus four of these coming from the edges, and four other reps $\delta^{j}_c$ (j runs from 1 to 4) belonging to the plaquettes themselves. Without any explicit computation, Witten now simply says that after summing the $\epsilon$ over all their indices (as required by the decomposition of a trace beforehand), we get a local factor associated to this vertex $G(\alpha_i,\delta^{(j)}_c,\epsilon)$, which is the 6j Wigner symbol (but he won't pause to show why).
I cannot find any source that would spell that relation out, i.e. show why we get precisely the 6j symbols in this computation (though their connection to the associator of the tensor product makes it plausible that we do). The real question is - the 6j symbol of what associator is this, and how would one go about and prove this?
I would be very grateful to anyone who can either explain this to me or direct me to a reference where this is discussed in more detail.
 A: Consider the finite dimensional unitary representations $\alpha,\beta,\gamma$ of the given compact group $G$ on corresponding vector spaces $V_1,V_2,V_3$. Let $|i\rangle_j,i=1,\dots,n_j$ be an orthnormal basis of $V_j$ where $dim V_j=n_j$. Then 
$\{|i\rangle_1\otimes|j\rangle_2\otimes|k\rangle_3\}$ forms an orthonormal basis of $V=V_1\otimes V_2 \otimes V_3$. An element $g\in G$ acts on the tensor product space $V$ as 
$$|i\rangle_1\otimes|j\rangle_2\otimes|k\rangle_3 \to \alpha(g) |i\rangle_1\otimes\beta(g)|j\rangle_2\otimes \gamma(g)|k\rangle_3 \tag 1$$ 
We can also find an orthogonal basis $e_{\mu},\mu=1,\dots,N$ (where $N=n_1n_2n_3$) of $V$ relative to which all elements $g \in G$ act as block diagonal matrices. More precisely, suppose with respect to basis $\{e_{\mu}\}$, the action of an element $g\in G$ on $V$ be denoted as 
$$e_{\mu}\to \rho (g)e_\mu \tag 2$$ 
then $\rho(g)^{\nu}_{\mu}=\langle \nu|\rho(g)|\mu\rangle$ is a block diagonal matrix where the dimensions of different blocks$^1$ are independent of $g$. Let 
$$|i\rangle_1\otimes|j\rangle_2\otimes|k\rangle_3=\sum_{\mu}\epsilon ^{\mu}_{ijk}e_{\mu}\tag 3$$
Acting with $g\in G$ on both sides of this equation gives 
$$\alpha (g)|i\rangle_1\otimes \beta (g)|j\rangle_2\otimes \gamma(g)|k\rangle_3=\sum_{\mu}\epsilon ^{\mu}_{ijk}\rho (g)e_{\mu} \tag 4$$
Taking the scalar product with $|i'\rangle_1\otimes |j'\rangle_2 \otimes |k'\rangle_3$ and using (3) we get 
$$\alpha(g)^{i'}_{i}\beta (g)^{j'}_{j}\gamma(g)^{k'}_{k}=\sum _{\mu,\nu} \rho^{\nu}_{\mu}(g)\epsilon^{*\nu}_{i'j'k'}{\epsilon}^{\mu}_{ijk}\tag 5$$
Now, according to the Peter-Weyl theorem (part 2), matrix elements of the irreducible representations of $G$ form an orthogonal basis of the space of square integrable functions on $G$ wrt the inner product 
$$(A,B)=\int_G dg\; A(g)^{*}B(g) \tag 6$$ 
where $dg$ is the Haar measure. So, if we integrate both sides of (5), the nonzero contribution on RHS will only come from the part of $\rho$ which is direct sum of identity representations. In other words, let $W\subseteq V$ be the subspace of $V$ on which $G$ acts trivially, and let $\{e_1,\dots e_m\}\subseteq $ $\{e_1,\dots e_m,\dots,e_N\}$ be the basis of $W$, then the integration of (5) gives 
$$\int_G dg\;\alpha(g)^{i'}_{i}\beta (g)^{j'}_{j}\gamma(g)^{k'}_{k}=\sum _{\mu,\nu =1}^{m} \delta^{\nu}_{\mu}\epsilon^{*\nu}_{i'j'k'}{\epsilon}^{\mu}_{ijk}=\sum _{\mu =1}^{m}\epsilon^{*\mu}_{i'j'k'}{\epsilon}^{\mu}_{ijk} \tag 7$$ 
where we have assumed that $Vol(G)=\displaystyle\int_G\; dg =1$ 
For the second part of your question, I would recommend these lecture notes. The  basic idea for computing Wilson loop averages is following -
For a surface with boundary, the partition function of two dimensional Yang-Mills theory depends on the holonomy along the boundary. Let the partition function of a surface of genus $h$ and one boundary be denoted as $Z_h(U,ag^2)$ where $U$ is the fixed holonomy along the given boundary, $a$ is the area of the surface and $g$ is the Yang-Mills coupling constant. Now, consider the simplest situation in which a Wilson loop $W$ in representation $R_W$ is inserted along a contractible loop $C$ on a closed surface of genus $h$, and area $a$. To compute the Wilson loop average, we first cut the surface along $C$, which gives a disc $D$ of area (say) $b$ and another surface $S$ of area $c=a-b$, genus $h$ and one boundary. Now the Wilson loop average is given by integrating over $G$ the product of i) the partition functions of $D$ ii) partition function of $S$ and iii) the trace of the Wilson loop in representation $R_W$ -
$$\langle W\rangle=\frac {1}{Z_h(ag^2)}\int \: dU \:Z_h(U,(a-b)g^2)\chi_{R_W}(U)Z_0(U^{-1},bg^2) \tag 8$$
The case of a self-intersecting Wilson loop too is not very different. 

$^1$ The smallest blocks form irreducible representations of $G$; Exactly which irreducible representations show up will depend on $\alpha,\beta,\gamma$; The same irreducible representations may also appear more than once
