system of homogeneous matrix equations Edit: Maybe I should make it clear that by a solution I mean a pair of matrices $(A,B)$ such that the identity below holds for any complex numbers $x,y$.
One of my friend asked me the following question (with some variation):
Let $A$,$B$ be two $n\times n$ complex matrices, let $x$ and $y$ be two complex variables. Suppose that 
$$
(xA+yB)^n=(x^n+y^n)I_n
$$
For all $x$ and $y$, Where $xA$ is the scalar multiplication of $A$ by $x$. $I_n$ is the identity matrix. 
Question: Can we get explicit solutions when $n$ is small? Can we at least say that the solutions forms a manifold of certain dimension for arbitrary $n$?
You can expand both sides, and move the right side to the left, getting a polynomial whose coefficients are matrices which is identical to $0$. This means the coefficients are all $0$, thus we get a system of homogeneous equations on $A$ and $B$. We can observe that the solution set is stable under taking conjugation by the same invertible matrix. Since $A^n=I_n$, we know that it's diagonalizable, whose eigenvalues are $n$-th roots of unity. Hence we can assume that $A$ is a diagonal matrix. By this simplification, using elementary methods, I was able to find the explicit solution for n=2 and 3. But as $n$ getting large, even for $n=5$, the redundant method requires one to solve a system of polynomial equations of 25 variables. I actually tried it using mathematica, but my poor computer became unresponsive very soon. Any suggestion is appreciated.
 A: Is there a reason why the power in the equation and the size of the matrices are the same? E.g., if you take two anti-commuting roots of $I_n$, you will always have $(xA+yB)^2 = (x^2+y^2)I_n$.
More generally, for two $n\times n$ matrices $A$ and $B$ with $A^k=B^k=I_n$ and $AB=qBA$ with $q$ a primitive $k$th root on unity, we have
$(xA+yB)^k = (x^k+y^k)I_n$.
This follows from the $q$-binomial theorem: if $AB=qBA$, then
$(A+B)^N = \sum_{\nu=0}^N C(N,\nu) A^\nu B^{N-\nu}$,
where the coefficients are determined by $C(N,0)=C(N,N)=1$ and the generalized "Pascal's triangle" recursion formula
$C(N+1,\nu) = C(N,\nu-1) + q^\nu C(N,\nu)$ for $1\le \nu \le N$. They are also called Gaussian binomial coefficients. 
A: A partial answer: Let $\Sigma$ denote the manifold $x^n+y^n=0$. Away from $\Sigma$, the equation and the fact that the roots of the polynomial $X^n-x^n-y^n$ are simple, tell you that $xA+yB$ is diagonalisable, with eigenvalues that differ from each other from a multiplicative $n$th root of unity. In addition, the multiplicity of an eigenvalue is locally constant.
The complement of $\Sigma$ in ${\mathbb C}^2$ is path-connected (true for every algebraic curve), but it is not simply connected. There are many non-trivial loops around $\Sigma$. When you follow an eigenvalue $\lambda(xA+yB)$ along a loop, you end with possibly different eigenvalue. Actually, there are so many loops that you find the following: if $\lambda$ is an eigenvalue of $xA+yB$, then $\omega\lambda$ is another one for every $n$th rooth of unity $\omega$. Eeven more: because the multiplicity remains constant when you follow the loop, $\lambda$ and $\omega\lambda$ have the same multiplicity.
In conclusion: for $(x,y)\not\in\Sigma$, $xA+yB$ is diagonalisable with eigenvalues at the vertices of a regular $n$-agon, and the eigenspaces have equal dimensions.
Since in your question you seem to take $n$ equal to the size of matrices, this means that the eigenvalues are simple away from $\Sigma$. Anyway, the existence of such a pair $(A,B)$ implies that $n$ (the power) divides $N$ (the matrix size).
