Yes. In general,
$$
\operatorname{vec}(A)^\top\operatorname{vec}(B) = \operatorname{tr}(A^\top B).
$$
It follows that
$$
v^\top (A \otimes B)v = \operatorname{vec}(V)^\top\operatorname{vec}(BVA^\top) = \operatorname{tr}(V^\top BVA^\top).
$$
By replacing $A$ by $\frac{1}{a}A$, we may assume that $a=1$ and $b>2$. With respect to the standard basis of $\mathbb R^{n\times n}$, the matrix $A$ represents the linear operator $f:\mathbb R^{n\times n}\to\mathbb R^{n\times n}$ defined by
$$
f(X)=\operatorname{tr}(X)I_n+(X+X^T)+(b-2)\operatorname{diag}(X).
$$
Hence $B$ represents the linear operator $g$ defined by
\begin{align}
g(X)&=Mf(MXM)M\\
&=M\left[\operatorname{tr}(MXM)I+MXM+(MXM)^T+(b-2)\operatorname{diag}(MXM)\right]M.\tag{0}
\end{align}
Observe that $g(X)=0$ when $X$ is skew-symmetric or when $MXM=0$. Hence
$$
K=\{X\in\mathbb R^{n\times n}:X=-X^T\}+\{X\in\mathbb R^{n\times n}:MXM=0\}\subseteq\ker g.
$$
As $g$ is self-adjoint (because $B$ is symmetric), in order to determine other eigenvectors or eigenvalues of $g$, it suffices to consider the restriction of $g$ on
$$
K^\perp=\{X\in\mathbb R^{n\times n}:X=MYM \text{ for some symmetric matrix } Y\}.
$$
Suppose $X\in K^\perp$. Then $MXM=M^2YM^2=MYM=X=X^T$. If this $X$ is an eigenvector $g$ corresponding to an eigenvalue $\lambda$, then
$$
\lambda X=g(X)
=\operatorname{tr}(X)M+2X+(b-2)M\operatorname{diag}(X)M.\tag{1}
$$
It follows that
\begin{align}
\lambda\|X\|_F^2
&=\langle X,\lambda X\rangle\\
&=\langle X,\,\operatorname{tr}(X)M+2X+(b-2)M\operatorname{diag}(X)M\rangle\\
&=\operatorname{tr}(X)\langle X,M\rangle+2\|X\|_F^2+(b-2)\langle X,M\operatorname{diag}(X)M\rangle\\
&=\operatorname{tr}(X)\operatorname{tr}(XM)+2\|X\|_F^2+(b-2)\operatorname{tr}(XM\operatorname{diag}(X)M)\\
&=\operatorname{tr}(X)\operatorname{tr}(XM^2)+2\|X\|_F^2+(b-2)\operatorname{tr}(XM\operatorname{diag}(X)M)\\
&=\operatorname{tr}(X)\operatorname{tr}(MXM)+2\|X\|_F^2+(b-2)\operatorname{tr}(MXM\operatorname{diag}(X))\\
&=\operatorname{tr}(X)\operatorname{tr}(X)+2\|X\|_F^2+(b-2)\operatorname{tr}(X\operatorname{diag}(X))\\
&=\operatorname{tr}(X)^2+2\|X\|_F^2+(b-2)\sum_{i=1}^nx_{ii}^2\\
&\ge2\|X\|_F^2.\tag{2}\\
\end{align}
Hence $g$ is positive definite on $K^\perp$, $\ker g=K$, every positive eigenvalue of $g$ is $\ge2$, and
$$
\operatorname{rank}(B)=\operatorname{rank}(g)=\dim K^\perp=\frac{1}{2}m(m+1),
$$
where $m$ denotes the rank of $M$. Moreover, $2$ is an eigenvalue of $g$ only if equality holds in $(2)$, i.e., only if $X$ is hollow ($\operatorname{diag}(X)=0$). This condition is also sufficient: if
$$
X\in S=\left\{X\in K^\perp:\operatorname{diag}(X)=0\right\}\tag{3}
$$
or equivalently, if $\operatorname{vec}(X)\in S_1=\operatorname{col}(M\otimes M)\cap\operatorname{col}(R)$ where $R$ is defined as in your question, then by $(0)$ we get $g(X)=2X$. Thus the multiplicity of the eigenvalue $2$ is equal to $\dim S\,(=\dim S_1)$.
This dimension can be zero. That is, $2$ is not necessarily an eigenvalue of $B$ or $g$. One obvious example is given by $M=0$, which makes $S=0$. For a non-trivial example, just try a numerical experiment that repeatedly samples some random orthogonal projectors $M\in\mathbb R^{3\times3}$ of rank $2$. Sometimes $2$ is an eigenvalue of $B$, but most of the time it isn't.
Since our problem formulation involves $\operatorname{diag}(X)$ or $\operatorname{diag}(MXM)$, it is not basis independent and we expect a nice answer when $M$ is, up to a permutation of rows and columns, in the form of $I_m\oplus0$ (where $m$ denotes the rank of $M$). Indeed, in this case one can easily see that the matrix subspace in $(3)$ consists of matrices of the form $F\oplus0$, where $F$ is any $m\times m$ hollow symmetric matrix. Therefore $2$ is an eigenvalue of multiplicity $\frac12m(m-1)$.
Best Answer
Not elegant but a direct proof. Fix some orthonormal $\{e_i\}$ and write $$ X = \sum_{ij} a_{ij} e_{i} e_{j}^*, $$ where $e_{i}e_{j}^*$ is the outer product of the basis elements $e_{i}$ and $e_{j}$. That is $e_{i}e_{j}^*$ is the matrix with a $1$ in the $(i,j)$ element and 0 elsewhere. Write the other matrices in a similar way then you can write $$ AXB^T = \sum_{ijkl} a_{ij} x_{jk} b_{lk} e_{i} e_{l}^*. $$ Note that the $\mathrm{vec}$ operation is linear and it acts on $e_{i}e_{j}^*$ as $\mathrm{vec}(e_{i}e_{j}^*) = e_{i} \otimes e_j$. Therefore, $$ \mathrm{vec}(AXB^T) = \sum_{ijkl} a_{ij} x_{jk} b_{lk} e_{i} \otimes e_{l}. $$
Then we can simply just calculate the RHS also by writing $$ A \otimes B = \sum_{ijkl} a_{ij} b_{kl} e_{i}e_{j}^* \otimes e_{k}e_{l}^* $$ and $$ \mathrm{vec}(X) = \sum_{mn} x_{mn} e_{m} \otimes e_n $$ which gives $$ \sum_{ijkl} a_{ij}b_{kl}x_{jl} e_{i} \otimes e_{k} $$ which is exactly what we had before up to relabelling $k \leftrightarrow l$.