Is there a formula or equation relating $\langle E\rangle$ and $\langle M\rangle$ (average spin per site) and $\langle E^2\rangle$ to temperature $T$ for the square lattice Ising model at zero magnetic field? I have done some simulations and know what the graphs are supposed to look like approximately, but is there some neat expression for them?
And is there known what the exact functions are for finite grid sizes?
Onsager computed the partition function of the 2D periodic square lattice (toroidal boundaries) Ising model. It is arguably one of the most elegant proof of modern statistical mechanics.
The original paper is available on the APS website below: (you will need institutional access)
L. Onsager, "Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition", Phys. Rev. (65), 1944. Link
Although I found it laying around on some university server: http://www.colorado.edu/physics/phys7230/phys7230_sp08/Onsager1944.pdf
In essence, he obtains the partition function $$Z(\beta,N,H=0)=(2\cosh(2\beta J) e^I)^N$$ with $$I=\frac{1}{2\pi}\int_0^\pi d\phi\ln\left(\frac{1}{2}\left[1+(1-\kappa^2\sin^2\phi)^{1/2}\right]\right)$$ where $$\kappa=\frac{2\sinh(2\beta J)}{\cosh^2(2\beta J)} $$
Traditionnal canonical ensemble techniques can be applied from there. Note that the free energy associated to $Z(\beta,N,0)$ is non analytic and that a phase transition arises when $\kappa=1$ . This correctly predicts that $T_c=\frac{2J}{k\ln(1+\sqrt{2})}$
Yes, these equations exist and can be derived from the partition function in JGab's answer.
The internal energy per spin is: $$u(\beta) = - \frac{\partial}{\partial \beta} \left( \ln(2) + \frac{1}{8 \pi^2} \int_{0}^{2\pi} \mathrm{d} q_1 \int_{0}^{2\pi} \mathrm{d} q_2 \\\ln \left[ \big( 1 - \sinh(2 \beta J) \big)^2 + \sinh(2 \beta J) \left( 2 - \cos q_1 - \cos q_2 \right) \right] \right)$$
The spontaneous magnetization per spin is: $$m_s(T) = \begin{cases} \left( 1 - \sinh^{-4} (2 \beta J) \right)^{\frac{1}{8}} & \text{for } T < T_c \\ \hskip 1.5cm 0 & \text{for } T \geq T_c \end{cases}$$ Note that this is strictly speaking not the magnetization $\langle m \rangle$ you asked for, but its limit $$m_s = \lim_{B \rightarrow 0^+} \lim_{N \rightarrow \infty} \langle m \rangle$$ where the two limits do not commute. See this question for a discussion of this important point.
