Equivalence classes of solutions to a matrix equation I am wondering whether the following problem can be mathematically formalized an solved. Namely, I want to quantify how many solutions there are to the equation
$$U_1 \cdot U_2 \cdot U_3 \cdot U_4 = U_4 \cdot U_3 \cdot U_2 \cdot U_1$$
where $\{U_a\}_{a=1}^4$ are unitary matrices of arbitrary (but the same) dimensions. For example, using the Pauli matrices, $U_1 = U_2 =\sigma_x$ with $U_3 = U_4 = \sigma_z $ would be a simple solution with $2\times 2$ matrices.
Clearly, there are infinitely many solutions, but one can organize them into enumerable equivalence classes. For example, replacing $U_a\mapsto U_a e^{\mathrm{i} \phi_a}$ maps one solution to another equivalent solution. Similarly, if $V$ is a unitary matrix of the same dimension, then $U_a \mapsto V\cdot U_a \cdot V^{-1}$ (with the same $V$ for all $a\in\{1,2,3,4\}$) gives another equivalent solution. However, not all solutions of the same dimension are equivalent; for example, the solution $U_1 = U_2 = \sigma_x$ and $U_3 = U_4 = \sigma_z$ cannot be related to the solution $U_1 = U_3 = \sigma_x$ and $U_2 = U_4 = \sigma_z$ using the listed equivalence transformations.
Finally, one can introduce the notion of "reducible" solutions, i.e. ones that are equivalent to $U_a$'s with a block-diagonal structure; and similarly, there are presumably "irreducible" higher-dimensional solutions. Clearly, the knowledge of irreducible solutions is sufficient for the identification of all possible solutions. The two example solutions listed above are both irreducible in this sense.
Thus, trying to be more precise now, I am wondering how to identify all equivalence classes of irreducible solutions to the equation at the beginning.
Remark: The fact that I chose to describe this problem using group-theoretical terminology is not an accident. The formulated problem should indeed be equivalent to looking for irreducible representations of the following free group with a constraint,
$$\Gamma = \left<a,b,c,d\,|\, a.b.c.d.a^{-1}.b^{-1}.c^{-1}.d^{-1}=1\right>.$$
This group corresponds to the maximal torsion-free subgroup of the Fuchsian symmetry group of the $\{8,8\}$-tessellation of the hyperbolic plane. Perhaps this offers some of you helpful geometric insights that I overlooked.
 A: The space you are actually interested in is known as the $U(n)$-character variety (there is a little lie here that I will ignore),
$$
X(\Sigma, U(n))= Hom(\pi_1(\Sigma), U(n))/U(n),
$$
where $\Sigma=\Sigma_2$ is the genus 2 closed oriented and connected surface and the quotient is taken with respect to the action of the group $G=U(n)$ on the space of homomorphisms via conjugation:
$$
g\cdot \rho(x)= g\rho(x) g^{-1}, x\in \pi_1(\Sigma).   
$$
Essentially all the discussion below only uses the fact that our surface has genus $p\ge 2$, so for now, that's all what I will be assuming.
By choosing $2p$ standard generators of $\pi_1(\Sigma)$,  the space $Hom(\pi_1(\Sigma), U(n))$ can be identified with a closed subset
of
$$
G^{2p}, G=U(n)
$$
(simply record the images $A_1, B_1,...,A_p,B_p$ of generators), subject to the equation
$$
[A_1,B_1]...[A_p,B_p]=1.
$$
There are two topologies one can consider in this setting, I will be always using the "classical" topology, coming from the Euclidean topology on the space of matrices. (The second natural topology is the Zariski topology.)
Since the space $Hom(\pi_1(\Sigma), U(1))$ is easy to understand (and you already do!) the study of the $U(n)$-character variety mostly reduces to that of the $SU(n)$-character variety.
The space of all homomorphisms (with values in $G=U(n)$ or $G=SU(n)$) is stratified as the union
$$
Hom^{irr}(\pi_1(\Sigma), G)\sqcup Hom^{red}(\pi_1(\Sigma), G).
$$
It is not hard to see that the former is open and, accordingly, the latter is closed. Similarly, we obtain a stratification of the character variety
$$
X=X^{red}\sqcup X^{irr}. 
$$
As it turns out, $Hom^{irr}(\pi_1(\Sigma), G)$ is a smooth submanifold of $G^{2p}$, the action of $G$ on it is free (and, of course, proper) and we obtain a principal $G$-fiber bundle
$$
Hom^{irr}(\pi_1(\Sigma), G)\to X^{irr}(\Sigma, G)
$$
and the latter is a smooth manifold.
Both facts (in greater generality) were first established by Andre Weil back in 1964 and were subsequently reproved  (invariably, in lesser generality) by various authors. Weil also gave the dimension formula which, in the special case at hand, reads:
$$
\dim X^{irr}(\Sigma, SU(n))= (2p-2) \dim(SU(n))= (2p-2)(n^2-1).
$$
Weil, André, Remarks on the cohomology of groups, Ann. Math. (2) 80, 149-157 (1964). ZBL0192.12802.
(Weil's writing is quite cryptic and, unless you already know what he is doing, is hard to follow.)
Weil did not know this, but $Hom(\pi_1(\Sigma), SU(n))$ has quadratic singularities at reducible representations and the subset of irreducible representations is dense in
$Hom(\pi_1(\Sigma), SU(n))$.
Personally, I like the proofs given by Goldman in
Goldman, William M., The symplectic nature of fundamental groups of surfaces, Adv. Math. 54, 200-225 (1984). ZBL0574.32032.
Goldman also proves the existence of a natural symplectic structure on the smooth manifold $X^{irr}(\pi_1(\Sigma), G)$.
Beyond these generalities, the study of $SU(n)$-character varieties is done using either the language of flat unitary connections, or of stable holomorphic vector bundles over surfaces, or of symplectic geometry.
The following notes by Liza Jeffrey:
Flat connections on 2-manifolds
while dated, serve as a good introduction. In particular, following her work with Jonathan Weitsman, she explains how to "parameterize" $X^{irr}(\Sigma, SU(2))$ using "Dehn coordinates" (based on a pair of pants decomposition of $\Sigma$) and then how to compute the cohomology rings of $SU(2)$-character varieties (this computation comes from formalization of some work of Ed Witten).
Jeffrey, L. C.; Weitsman, J., Toric structures on the moduli space of flat connections on a Riemann surface: Volumes and the moment map, Adv. Math. 106, No. 2, 151-168 (1994). ZBL0836.58004.
Jeffrey, Lisa C.; Weitsman, Jonathan, Bohr-Sommerfeld orbits in the moduli space of flat connections and the Verlinde dimension formula, Commun. Math. Phys. 150, No. 3, 593-630 (1992). ZBL0787.53068.
See also this paper:
Charles, Laurent; Marché, Julien, Multicurves and regular functions on the representation variety of a surface in (SU(2)), Comment. Math. Helv. 87, No. 2, 409-431 (2012). ZBL1246.57022.
(a free pdf file can be found here).
Such a parameterization is the closest thing available to your request for a "description of all irreducible representations."
In the case of the genus 2 surface you are interested in, there are two pair of pants decompositions of $\Sigma$; both are along three simple loops $a, b, c$ in $\Sigma$. Each equivalence class $[\rho]\in X^{irr}(\Sigma, SU(2))$ is uniquely determined by two sets of real numbers assigned to $a, b, c$: The first set of numbers is derived from the traces $tr(\rho(a)), tr(\rho(b)), tr(\rho(c))$ of the matrices in $SU(2)$ and they satisfy a set of "simplex inequalities" (which can be interpreted as "triangle inequalities" for spherical triangles). The second set of numbers are "twist" parameters; they are unconstrained. (There are an extra complications here coming from the fact that an irreducible representation can become reducible when restricted to a pair of pants.)
Similar parameterizations exist for $n\ge 3$ but everything gets much harder (as $n$ increases). The "simplex inequalities" defining the moment polytope were generalized by Belkale  to a much more complex set of linear inequalities:
Belkale, Prakash, Local systems on ({\mathbb{P}}^1-S) for (S) a finite set, Compos. Math. 129, No. 1, 67-86 (2001). ZBL1042.14031.
These inequalities are in terms of "quantum Schubert calculus," but do admit an inductive (in terms ofn $n$) description also given by Belkale.
I'll stop here since the answer is already quite long.
