In reply to Rick Wicklin.

I am not totally sure, other than to say that is to do with the pattern of eigenvalues. I read "Odd Magic Powers", A. C. Thompson, The American Mathematical Monthly, 101(1994), pp. 339-342, and then started experimenting in IML. Only odd powers will work and it seems that squares of order 4n produced by your function are specially symmetrical (but not pan-diagonal).
.
I have an even more fantastic example. This produces a set of 6 magic squares!
.
x = magic(8);
print x [f=3.0];
do cycle=1 to 6;
x = rank( mod( x**5, 346) );
print x [f=3.0];
end;
.
If you try larger squares or larger powers, be careful not to exceed the largest integer that can be exactly represented.