In the TQFT language, quasiparticles correspond to Wilson loop operators. It is well-known that quasiparticles can have non-trivial braiding statistics.
Take $2+1$ dimensional Abelian Chern-Simons theory with $U(1)$ gauge group as an example,
\begin{equation}
S_{CS} = \int \frac{k}{4\pi}A\wedge dA - q A\wedge *j
\end{equation}
where $j$ is a source term representing quasiparticle currents (which carries $q$ units of $A-$gauge charge).
The equation of motion reads
\begin{equation}
\frac{2\pi q}{k}*j = dA
\end{equation}
This implies that each quasiparticle not only carries $q$ units of $A-$gauge charge, but are dressed with $\frac{2\pi q}{k}$ units of $A-$gauge flux. So braid a charge$-q$ particle around a charge$-q^{\prime}$ particle induces (by the Aharonov-Bohm effect) a non-trivial $U(1)-$phase $\theta = \frac{2\pi}{k}q q^{\prime}$.
Wilson lines in TQFT captures the above statistics nicely. In order to see this, one can calculate the regulated expectation value of two line operators,
\begin{equation}
<\exp(iq\oint_{C_1}A)\exp(iq^{\prime}\oint_{C_2}A)> = <\exp(\frac{2\pi i}{k}q q^{\prime}L(C_1,C_2))>
\end{equation}
where $L(C_1,C_2)$ is the linking number of the two lines. The process of braiding one charge around another corresponds to $L=1$. The Chern-Simons theory thus automatically calculates the linking number of chosen line operators and multiplies the wavefunction by the appropriate phase.
Another example is $2+1$ dimensional BF theory,
\begin{equation}
S_{BF} = \int \frac{k}{2\pi}B\wedge dA - A\wedge *j - B\wedge *J
\end{equation}
where both $A$ and $B$ are $1-$form gauge fields, and $j$ and $J$ are quasiparticle current carrying $A$ and $B$-gauge charges respectively.
The equation of motion reads
\begin{align}
dA &= \frac{2\pi}{k}*J \nonumber \\
dB &= \frac{2\pi}{k}*j
\end{align}
The magnetic flux of $A$ is attached to the charges of $B$ and vice versa, and braiding $q$ units of $A-$charge around $m$ units of $B-$charge induces a phase $\frac{2\pi q m}{k}$. This is again nicely captured by the following equation
\begin{equation}
<\exp(iq\oint_{C_1}A)\exp(im\oint_{C_2}B)> = <\exp(\frac{2\pi i}{k}qmL(C_1,C_2))>
\end{equation}
My question is, are there any similar characterizations for braiding statistics in $3+1$ dimension? In $3+1$ dimension, generic excitations are quasiparticles as well as quasistrings. So we need to be able to characterize braidings between Wilson loops and Wilson surfaces. I think equation (2.13) on page 8 of http://arxiv.org/abs/1011.5120 defines something similar.
However, besides braidings between quasiparticles and quasistrings. The really fundamental braiding process in $3+1$ dimensional gauge theories are the so-called three-loop braiding process, as argued in http://arxiv.org/abs/1403.7437. The statistics is between 3 loops carrying gauge fluxes, where two loops braid with each other while both are linked to a third base loop. I'm wondering if there is an analogous characterization of the above process in terms of Wilson surfaces. Is there something called "linking number" of Wilson surfaces? If yes, how is it defined?