The $SU(3)_1$ WZW model possesses three primary operators, namely $\mathbf{1}$, $\mathbf{3}$ and $\bar{\mathbf{3}}$ which are identified as the identity matrix, $g(z,\bar{z})$ and $g^{-1}(z,\bar{z})$ with $g \in SU(3)$ . Now, the model was famously introduced in the context of non-abelian bosonization and thus has a description in terms of 3 free fermions with gauged $U(1)$ symmetry, see e.g. 1 and 2 for the construction. The upshot is that (staying in the notation of 1), using the free fermions (sum over $i=1,2,3$ suppressed)
\begin{align}
S[\chi, \chi^\dagger]=\frac{1}{\pi} \int\left(\chi_2^{\dagger} \bar{\partial} \chi_1+\chi_1^{\dagger} \partial \chi_2\right) d^2 x,
\end{align}
we can identify fermionic bilinears as
$$
g^{i j}(z, \bar{z}) \sim \chi_2^i(\bar{z}) \chi_2^{j \dagger}(z) e^{-\varphi(z,\bar{z})}, \qquad (g^\dagger)^{i j}(z, \bar{z}) \sim \chi_1^i (z) \chi_1^{j \dagger} (\bar{z}) e^{\varphi(z,\bar{z})}
$$
where normal ordering is implied and the normalization of the action of the free boson $\varphi(z,\bar{z})$ is chosen such that the vertex operators $e^{\pm\varphi(z,\bar{z})}$ have conformal dimension $(h, \bar{h})=(-1/6,-1/6)$.
My goal is now to calculate the $gg$ OPE in Fermionic language and to somehow recover the $\mathbf{3} \times \mathbf{3} \to \bar{\mathbf{3}}$ fusion rules. This means that we need to recover (focusing on the holomorphic sector) $$ \chi_2^{i \dagger}(z) \chi_2^{j \dagger}(w)e^{-\varphi(z)}e^{-\varphi(w)} \sim \frac{T_{ijk}}{(z-w)^{1/3}} \chi_1^{k}(w)e^{\varphi(w)} + (\text{non-singular}), $$ where $T_{ijk}$ is some tensor to be determined. Of course one also needs to include the anti-holomorphic sector for this identity to make sense (otherwise we would change Grassmann parity).
My question would now be: How to derive this OPE?
Some things I figured out: Since the LHS of the OPE is anti-symmetric in $i$ and $j$, the RHS must is an anti-symmetric tensor, meaning it can be dualized to $\epsilon_{ijk}O^k(w)$ with $O^k(w)$ some three-component operator. Now how should $O^k$ transform under $SU(3)$? A nice identity worked out here is $\epsilon_{ijk} (Ua)^i (Ub)^j c^k=\epsilon_{ijk} a^i b^j (U^\dagger c)^k$ for $U \in SU(3)$. Applied to our problem, this tells us that $O$ transforms with $U^\dagger$, which looking at $S[\chi,\chi^\dagger]$ implies that it must be $O^k=\chi_1^k$. Thus, we have found our operator.
To make it more concrete, in eq. (30) of 1 it is claimed that (they use the anti-hol sector) $$ \chi_2^{i \dagger}(z) \chi_2^{j \dagger}(w)e^{-\varphi(z)}e^{-\varphi(w)} \sim {(z-w)^{-1/3}} \underbrace{:\chi_2^{i \dagger}(z) \chi_2^{j \dagger}(w):e^{-2 \varphi(w)}}_{= g^{-1}_{ij}} + (\text{non-singular}) $$ where the singularity comes from the vertex operator OPE. Thus, we should be able to identify $$ \epsilon_{ijk} \,\chi_1^{k}(w)e^{\varphi(w)}=:\chi_2^{i \dagger}(z) \chi_2^{j \dagger},(w):e^{-2 \varphi(w)} $$ including the OPE-ceofficient $C_{\mathbf{3}\mathbf{3}\mathbf{\bar{3}}}=1$.
How can I see this?
Many thanks in advance!
