I asked a similar question here When does $AA^T$ commute with $A^T$?.
Let $A \in M_n(\mathbb R)$ with spectral radius $\rho(A) < 1$. So the series $X = \sum_{j=0}^{\infty} A^j(A^T)^j$ makes sense. I want to determine when $XA^T = A^T X$.
If $A$ is normal, this will certainly hold since $A^j (A^T)^j A^T = A^T A^j (A^T)^j$ for all $j \in \mathbb N$ as $A^T$ is polynomial in $A$. But I am wondering whether there is a larger class of $A$ resulting this commutativity.
If $X$ and $A^T$ commute, $A$ must be normal.
Note that when $\rho(A)<1$, the infinite series always converges to the unique symmetric positive definite solution to the Lyapunov equation $AXA^T=X-I$ (see here for the reason). Since $X-I$ is positive semidefinite, by a change of orthonormal basis, we may assume that $X=\lambda_1I_{k_1}\oplus\cdots\oplus\lambda_mI_{k_m}$ for some distinct real numbers $\lambda_1,\ldots,\lambda_m\ge1$. So, if $A^T$ commutes with $X$, it must be a block diagonal matrix of the form $Q_1\oplus\cdots\oplus Q_m$, where each $Q_j$ has the same size as $I_{k_j}$. But then the equality $AXA^T=X-I$ implies that $Q_j^TQ_j=\frac{\lambda_j-1}{\lambda_j}I_{m_j}$. Hence all $Q_j$s and in turn $A$ are normal.