Number of solutions $X$ to $AX=XB$ in $\mathbb F_2$ It is a well-known theorem that in an arbitrary field $F$, if $A$ is an $m\times m$ square matrix and $B$ is an $n\times n$ square matrix, then there is a unique $m\times n$ solution $X$ to the equation 
$$AX=XB$$
if and only if $A$ and $B$ share no eigenvalues.
My question is this: what can be said about the number of solutions $X\in \mathbb F_2$ to the above equation? Is the number of nonzero solutions equal to the number of common eigenvalues? If so, how can one prove this?
 A: $X$ maps each eigenvector of $B$ to either $0$ or an eigenvector of the same eigenvalue of $A$, since if $Bv = \lambda v$, $AX v = X B v = \lambda X v$.
Let's suppose for simplicity that $B$ is diagonalizable, with a basis $b_j$ of eigenvectors corresponding to eigenvalues $\lambda_j$ that may or may not be eigenvalues of $A$.
For each $j$, we can choose $X b_j$ to be any member of the eigenspace of $A$ for $\lambda_j$, and once we have done that we have specified a solution $X$.  Thus the number of solutions for $X$ is the product over $j$ of the cardinalities of the eigenspaces of $A$ for $\lambda_j$ (note that this is $1$ if $\lambda_j$ is not an eigenvalue of $A$).
A: Your statement of well-known theorem requires saying that the eigenvalues are in the algebraic closure of the given field $F$. This is an easy version of Cecioni-Frobenius theorem suggested in the comments by 'user'. 
Method 1 : Module Theory
Under your assumption for matrices $A$, $B$, $X$, define 
$$
C_{A,B}=\{X\in M_{m\times n}(F) \mid AX - XB = 0 \}.$$
Let $\nu_{A,B}$ be the dimension of $C_{A,B}$ over $F$. The actual version of Cecioni-Frobenius theorem gives a formula
$$
\nu_{A,B}=\sum_p(\deg p)\sum_{i,j} \min\{\lambda_{p,i}, \mu_{p,j}\}
$$
Here the first sum is over all irreducible polynomials $p$ which are common in the primary decompositions of $F[x]$-modules $M^A$ and $M^B$, and the indices $i, j$ of second double sum is from the partition $\lambda_p=\sum_i \lambda_{p,i}$ that indicates the powers of $p$ in $p$-primary part of $M^A$, $\mu_p=\sum_j \mu_{p,j}$ that of powers of $p$ in $p$-primary part of $M^B$.
Then the number of solutions to $AX=XB$ is $2^{\nu_{A,B}}$. 
Method 2 : Kronecker Product
Alternatively, there is a module-theory-free way. You will have an advantage of not having to solve for eigenvalues of $A$, $B$. But, a disadvantage is that there is almost no hope for writing down an explicit formula for $\nu_{A,B}$. This is by Kronecker product, the method is mainly this: Solution of a Sylvester equation?
In our case, the equation $AX-XB=0$ can be written as 
$$
(I_n \otimes A - B^T \otimes I_m) \textrm{vec}(X)=0
$$
where $\textrm{vec}(X)$ is a stack of columns of $X$ in order. 
Then perform elementary row operations on the $mn\times mn$ matrix  $I_n \otimes A - B^T \otimes I_m$. Find the number of free-variables, which is $\nu_{A,B}$. Again the number of solutions to $AX=XB$ is $2^{\nu_{A,B}}$. 
A: To ask for the number of solutions of $f(X)=AX-XB=0$ is equivalent to ask for $dim(\ker(f))$. Moreover, the required dimension is the same over the underlying field $K$ or over its algebraic closure $\overline{K}$. 
If we stack the considered matrices row by row, then $f=A\otimes I-I\otimes B^T$. Let $(\lambda_i),(\mu_i)$ be the spectra of $A,B$ in $\overline{K}$. Since the functions $X\mapsto AX,X\mapsto XB$ commute, $spectrum(f)=(\lambda_i-\mu_j)_{i,j}$. Unfortunately, that does not give the required dimension; we can only say that $\ker(f)\not= 0$ iff $degree(gcd(\chi_A,\chi_B))\geq 1$, where $\chi_.$ is the characteristic polynomial.
The Cecioni-Frobenius (CF) theorem (1910) provides a solution if we can calculate the invariant factors of $A,B$.
$dim(\ker(f))=\nu_{A,B}=\sum_{i,j}degree(gcd(d_i(A),d_j(B))$ where $d_i$ is the $i^{th}$ invariant factor. 
I had heard about this theorem in 2009 and then I had completely forgotten about it. It was used by Byrnes-Gauger (1979) to show that
$A,B$ are similar iff $2\nu_{A,B}=\nu_{A,A}+\nu_{B,B}$.
Note that, in great dimension, the invariant factor algorithm (the calculation of the Smith normal form of $xI-A$) is not effective. Fortunately, Byrnes-Gauger before proved (1977) this extraordinary result
$A,B$ are similar iff $\det(xI-A)=\det(xI-B),rank(A\otimes I-I\otimes A)=rank(B\otimes I-I\otimes B)=rank(A\otimes I-I\otimes B)$.
Here is an example of using the CF theorem
let $J_k$ be the nilpotent Jordan block of dimension $k$, $A=diag(J_3,J_2),B=diag(J_4,J_1)$. Then the Smith forms are 
$S(A)=diag(1,1,1,x^2,x^3),S(B)=diag(1,1,1,x,x^4)$.
Then $\nu_{A,B}=deg(gcd(x^2,x))+deg(gcd(x^2,x^4))+deg(gcd(x^3,x))+deg(gcd(x^3,x^4))=7$.
A: This is just a special case of Sylvester equation.
Let $A$ and $B$ be operators with spectra $\sigma(A)$ and $\sigma(B)$, respectively.  If $\sigma(A)$ and $\sigma(B)$ are disjoint, then the equation
$AX-XB=Y$ has a unique solution $X$ for every $Y$.  Moreover, the solution is given by $X=\sum_{n=0}^\infty A^{-n-1}YB$.  Hence, for $Y=0$ unique solution is $X=0$ as Robert Israel claimed.
The proof of the theorem above can be found in Bhatia's Matrix Analysis.
A: Suppose $X_1$ and $X_2$ are both solutions. Then $A(X_1+X_2)=AX_1+AX_2=X_1B+X_2B=(X_1+X_2)B$. Thus, $X_1+X_2$ is a solution. So the number of solutions can be $0$, $1$, or infinity. This is a common situation in linear algebra: you either have no solutions, a unique solution, or infinite solutions.
