I think the problem is that M and M.M both have the eigenvalue 1 with multiplicity 2 or higher (the multiplicity of 1 for M is 2 while it is 3 for M.M).

That means that the eigenvectors to be returned by Eigensystem belonging to eigenvalue 1 are not uniquely defined - any orthogonal basis of the eigenspace of eigenvalue 1 would do.

I think the problem is that M and M.M both have the eigenvalue 1 with multiplicity 2 or higher (the multiplicity of 1 for M is 2 while it is 3 for M.M).

That means that the eigenvectors to be returned by Eigensystem belonging to eigenvalue 1 are not uniquely defined - any orthogonal basis of the eigenspace of eigenvalue 1 would do.

Moreover, the eigenvalues of M.M are sorted in a different order than those of M:

M = {{0, 1, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 
    1}, {0, 0, 1, 0, 0}};
{λ, P} = Eigensystem[M];
P = Transpose[P];

{λp, Pp} = Eigensystem[M.M];
Pp = Transpose[Pp];
Max[Abs[λ^2 - λp[[{1, 2, 3, 5, 4}]] // N]]

1.11022*10^-16