Let us analyze the case where M is diagonalizable, since the defective case of a non-diagonalizable matrix (Jordan blocks) can be obtained at the end as the limit of diagonalizable matrices.

The diagonalization is M = R\Lambda R^{-1} where R is a matrix whose columns are (right) eigenvectors. Then M^T = (R^{-1})^T \Lambda R^T, but this is also a diagonalization and so we see by inspection that (R^{-1})^T = L, a matrix whose columns are *left* eigenvectors. (This is one proof of the general property that we can choose L^T R = I, i.e. left and right eigenvectors can be chosen orthonormal under the “unconjugated inner product” x^T y.)

That is, M = R \Lambda L^T with L^T R = I, and it follows that M^n = R \Lambda^n L^T.

In your case, all of the eigenvalues have |\lambda|<1 except for eigenvalues with \lambda=1, so the limit \Lambda^\infty is a diagonal matrix that projects onto the rows/columns corresponding to the \lambda=1 eigenvalues. Hence M^\infty = R \Lambda^\infty L^T = R_0 L_0^T, where L_0 and R_0 are the left and right eigenvectors for \lambda=1, equivalent to left and right nullspaces \hat{L}_0 and \hat{R}_0 of M-I. Technically, the left nullspace is the nullspace of M^* - I, but this conjugation is cancelled if we also use M^\infty = \hat{R}_0 \hat{L}_0^* instead of the transpose.

However, we have to be careful, especially in the case where \lambda=1 has multiplicity > 1, since there are multiple possible choices of the left and right eigenvectors/nullspaces and we have to choose them consistently to ensure that \hat{L}_0^* \hat{R}_0 = I. If we choose them arbitrarily, we might in general need some change of basis \hat{R}_0 \to \hat{R}_0 B, giving M^\infty = \hat{R}_0 B \hat{L}_0^* for some (invertible) B. We can easily determine B by the fact that M^\infty \hat{R}_0 = \hat{R}_0, which immediately gives B = (\hat{L}_0^* \hat{R}_0)^{-1}.

So, in the end, for any diagonalizable matrix M with the spectrum you assumed, we have M^\infty = \hat{R}_0 (\hat{L}_0^* \hat{R}_0)^{-1} \hat{L}_0^*, in terms of left and right nullspace bases \hat{L}_0 and \hat{R}_0 of M-I (regardless of how they are obtained/normalized).

For a defective matrix M, as long as the \lambda=1 eigenvalue is not itself defective, if we take the defective matrix M as the limit of diagonalizable matrices (which is always possible since diagonalizable matrices are a dense subset of all square matrices), we get the same answer (the limit of the \lambda=1 eigenvectors still span the nullspaces of M-I).

PS. I guess you get this spectrum from some kind of Markov property?