Braiding statistics of anyons from a Non-Abelian Chern-Simon theory Given a 2+1D Abelian K matrix Chern-Simon theory (with multiplet of internal gauge field $a_I$) partition function:
$$
Z=\exp\left[i\int\big( \frac{1}{4\pi} K_{IJ} a_I \wedge d a_J + a \wedge * j(\ell_m)+ a \wedge * j(\ell_n)\big)\right]
$$
with anyons (Wilson lines) of $j(\ell_m)$ and $j(\ell_n)$.
One can integrate out internal gauge field $a$ to get a Hopf term, which we interpret as the braiding statistics angle, i.e. the phase gained of the full wave function of the system when we do the full braiding between two anyons:
$$
\exp\left[i\theta_{ab}\right]\equiv\exp\left[i 2 \pi\ell_{a,I}^{} K^{-1}_{IJ} \ell_{b,J}^{}\right]
$$
see also this paper and this paper.
I would like to know the way(s) to obtain braiding statistics of anyons from a Non-Abelian Chern-Simon theory? (generically, it should be a matrix.) How to obtain this braiding matrix from Non-Abelian Chern-Simon theory?
 A: How to obtain this braiding matrix from Non-Abelian Chern-Simon theory?
To obtain braiding matrix $U^{ab}$ for particle $a$ and $b$, we first need to know the dimension of the matrix. However, the dimension of the matrix for Non-Abelian Chern-Simon theory is NOT determined by $a$ and $b$ alone. Say if we put four particles $a,b,c,d$ on a sphere, the dimension of the degenerate ground states depend on $a,b,c,d$. So even the dimension of the braiding matrix $U^{ab}$ depends on $c$ and $d$. The "braiding matrix" $U^{ab}$ is mot deterimened by the two particles $a$ and $b$.
Bottom line: physically, the Non-Abelian statistics is not described by the "braiding matrix" of the two particles $a$ and $b$, but by modular tensor category.
A: The (unitary) "phase" factor for non-Abelian anyons satisfies the
(non-Abelian) Knizhnik-Zamolodchikov equation:
$$\big (\frac{\partial}{\partial z_{\alpha}} + \frac{1}{2\pi k} \sum_{\beta \neq \alpha} \frac{Q^a_{\alpha}Q^a_{\beta}}{z_{\alpha} - z_{\beta}}\big )U(z_1, ....,z_N)  = 0 $$
Where $z_{\alpha}$ is the complex plane coordinate of the particle
 $\alpha$ , and $Q^a_{\alpha}$ is the matrix representative of the $a-$th
 gauge group generator of the particle $\alpha$ and $k$ is the level .
Please, see the following two articles by Lee and Oh (article-1,
article-2).
In the first article they explicitly write the solution in the case of
the two-body problem:
$$U(z_1, z_2)  = exp( i\frac{Q^a_1Q^a_2}{2\pi k} ln(z_1-z_2))$$
The articles describe the method of solution:
The non-Abelian phase factor can be obtained from a quantum mechanical
model of $N$ particles on the plane each belonging possibly to a
different representation of the gauge group minimally coupled to a gauge field with a Chern-Simons term in the Lagrangian. 
The classical field equations of the gauge potential can be exactly solved and
substituted in the Hamiltonian. The reduced Hamiltonian can also be exactly solved. Its solution is given by the action of a unitary phase factor on a symmetric wave function. This factor satisfies the Knizhnik-Zamolodchikov equation. 
The unitary phase factor lives in the tensor product Hilbert space of
the individual particle representations. The wave function is a vector in this Hilbert space valued holomorphic function depending on the $N$ points in the plane.
