H[α_][M_]:=-Sqrt[Length[M]]Total[Flatten@((Abs[M])^α)]

metroStep[β_,H_,δ_][M_]:=If[RandomReal[]<Exp[-β(H[#]-H[M])],#,M]&@(RotationMatrix[δ,{RandomVariate[NormalDistribution[],Length[M]],RandomVariate[NormalDistribution[0,δ],Length[M]]}].M)

metroRun[n_][β_,H_,δ_][M_]:=Nest[metroStep[β,H,δ],M,n]

data={metroStep[1,H[1],0.2][IdentityMatrix[4]]}

data=NestList[metroRun[40][10,H[2],0.5],data[[-1]],5000][[5;;]];

Histogram[H[2]/@data/(4^2)]

1+Mean[H[1]/@data]/(4^2)

H[1][IdentityMatrix[12]]/(12^2)//N

t1=FullSimplify[Total[Flatten[((Abs[RotationMatrix[θ]])^(2-ϵ))]],{θ∈Reals,ϵ∈Reals}]

FullSimplify[Series[t1,{ϵ,0,2}],{θ∈Reals}]

T[n_]:=Total[ArrayPad[RotationMatrix[θ[#]],{{2(#-1),2(n-#-1)},{2(#-1),2(n-#-1)}}]&/@Range[n]]

((Abs[(Transpose[Array[p,{10,10}]].T[6].Array[p,{10,10}])])^(2-ϵ))//Total