The configuration space $\mathcal {M}$ of $N$ non-Abelian anyons on the $2-$ sphere $S^2$ consists of $N$ copies of $S^2$ excluding the collision points and in addition $N$ irreducible representations of a non-Abelian Lie group $G$ corresponding to the highest weight vectors $j_k$, $k= 1, .,.,.,N$.
The corresponding (physical) Hilbert space is given by:
$$H = \mathrm{inv}\left(\bigotimes_{k=1}^{N} H^{(j_k)}\right)$$
Where, $H^{(j_k)}$ is the Hilbert space of the $k-$th anyon carrying a representation $j_k$ of a non-Abelian Lie group $G$. ($j_k$ is the highest weight vector of the representation).
The physical Hilbert space includes only the scalar representations in the tensor product .
The braiding matrix is a non-Abelian holonomy $B$ of a flat connection (which will be described below) on the trivial bundle $\mathcal {P} = \mathcal {M} \otimes G$, corresponding to a closed path in the configuration space in which one anyon say the $m-$th one encircles a single other anyon say the $n-$th. (The exact trajectory is immaterial since the connection is flat). In this case the braiding matrix $B = q^t$ is given by: (Here I am writing explicitly the short hand notation of Pachos):
$$t = \sum_{a=1}^{\dim(G)} I_1 \otimes . . . \otimes T^a_{j_m}\otimes . . . \otimes T^a_{j_n}\otimes ... \otimes I_k$$
The sum in the braiding matrix is over an orthonormal basis of the Lie algebra $\mathfrak{g}$ of $G$. $I_k$ are unit matrices, so the braiding matrix acts only on the tensor product of the local Hilbert spaces corresponding to the anyons $m$ and $n$.
Please notice that due to fact that $t$ is invariant under $G$, the braiding works within the physical Hilbert space of invariants.
The braiding matrix can be obtained from a quantum mechanical model of the Anyon system. I am following here Oh and Verlinde . A system of particles in $2+1$ dimensions acquires anyon statistics when it is coupled to a Chern-Simons term. In the case of non-Abelian Anyons, the particles need to have internal structure and to be coupled to a non-Abelian Chern Simons term. Oh works with the following theory of non-relativistic Anyons whose action (for $SU(2)$ anyons is given by:
$$ I = \int dt \sum_{k=1}^N \left( \frac{m_k}{2} \dot{q_k}^2 - Q_k^a (\theta_k, \phi_k)\left(A_i^a (q_k, t)\dot{ q_k }^i - A_0^a(q_k, t)\right) +j_k\cos(\theta_k)\dot{\phi_k}\right ) + \quad k \int dt d^2x \epsilon^{\mu \nu \lambda} (A_{\mu }\partial_{\nu} A_{\lambda } + \frac{2}{3} A_{\mu }A_{\nu} A_{\lambda } )$$
This is equation (2) in the article with a slight change of notation. $j_k$ is the spin of the $k-$th anyon; the second term in the action is the non-Abelian Lorentz term interaction. The third term is the symplectic potential of the two-spheres consisting of the classical phase spaces of the spin degrees of freedom. $ Q_k^a$ are the classical non-Abelian charges, which due to the fact that the group is $SU(2)$ are given by the components of the angular momentum (The Poisson algebra of the sphere):
$$(Q_k^1, Q_k^2, Q_k^3) = j_k \left(\cos(\phi_k) \sin(\theta_k), \sin(\phi_k) \sin(\theta_k), \cos(\theta_k) \right)$$
The Gauss' law (equations of motion of $A_0^a$) takes the form:
$$k F_{12}^a = \sum_{k=1}^N Q_k^a \delta^2(x-x_k)$$
Except at the particle points the gauge field is flat, and the flat connection solving the Gauss Law is given in the holomorphic gauge ($z = x^1+ix^2$)
$$A_{z}^a = \sum_{k=1}^N Q_k^a \frac{1}{z-z_k}$$
The charges $ Q_k^a $ quantize into the representation matrices $T_k^a$ of the $\mathfrak{su}(2)$ algebra. The holonomy of the connection becomes:
$$B = Pe^{\frac{1}{2 \pi k}\int_{\Gamma} \sum_k dz^k \sum_{l\ne k} dz^l T_k^a \otimes T_l^a \frac{1}{z_k-z_l}}$$
In the two-body case the exponent becomes self commuting over time and the path ordering becomes unnecessary and we get the final result:
$$B = e^{\frac{i}{2k} T_1^a \otimes T_2^a}$$
