La respuesta a esta OP está en el artículo citado por @loup blanc.
He añadido este post sólo para completar el inicial de la pregunta con una respuesta que sea fácilmente accesible a los lectores, que resume los principales resultados en el papel.
Parece que la afirmación de que una ecuación de la forma:
$$
(1) \qquad A_nX^n+A_{n-1}X^{n-1}+\cdots + A_0=0
$$
(donde $A_i$ $X$ son cuadradas complejas matrices de orden $k$) $\binom {nk}{k}$ soluciones, fue encontrada por J. J. Sylvester, pero fue comunicada sin una prueba.
La prueba se administra en el citado trabajo se demuestra que las matrices $X_h$ que son soluciones de la ecuación de $(1)$ son matrices que tienen como valores propios $k$ números de $\{\lambda_i\}\quad (\lambda_1 \in \{1,\cdots,nk\})$ , donde el $\lambda_i$ son soluciones de la ecuación:
$$
(2) \qquad \mbox{det}\left(\lambda^n A_n+ \cdots+\lambda A_1+A_0\right)=0
$$
que es, obviamente, una ecuación de grado $nk$, por lo que el $X_h$ son exactamente $\binom {nk}{k}$.
En el documento también se demostró que el grupo de Galois de $(2)$$S_{nk}$, por lo que las entradas de la solución de $(1)$ puede ser calculado por radicales si y sólo si $n=k=2$.