# Estimating the heat capacity of ising model

I am have written a Metropolis-Hastings algorithm and am currently trying to compare it to the analytical results for the 2D Ising model. The free energy seems reasonable but the heat capacity I'm getting is way off.

I am using the equation $$C_{v}=\frac{\beta}{T}[\left-\left^{2}]$$. Using $$\beta=1$$ would imply that $$T$$ is on the order of $$10^{23}$$ so the prefactor becomes really small. My understanding was that the variance in the energy will be humongous at such high temperatures but my variance is not that large. I think it is because I am using a $$16\times16$$ lattice and the variance on a finite system (I'm guessing) is bounded. But the computation time is already really slow (5 minutes for $$10^{6}$$ iterations of Monte Carlo).

I have tried using $$C_{v}=-\beta^{2}\frac{\partial^{2} f}{\partial \beta^{2}}$$ with a numerical approximation for the second derivative. But my free energy values are also numerically calculated so I sometimes get nonsensical results like negative heat capacity. I was wondering if there are other ways of estimating the heat capacity. Or if I'm stuck with using the variance of the energy, how large a lattice do I need to have?

Above the critical temperature, the variance of the energy diverges like the number $$N$$ of spins in your system. It would be far more relevant to estimate the specific heat $$C_v/N$$, for which one has an exact expression in the thermodynamic limit (see Chapter 5 of McCoy and Wu's book or this page).
More precisely, I would suggest that you measure directly the variance of the energy or, better, of the energy density $$H/N$$ (where $$H$$ is the Hamiltonian). The latter can then be compared with the exact solution since $$\operatorname{Var}_\beta(H/N) = \frac{d^2}{d\beta^2} \frac{1}{N}\log Z_{N,\beta}$$, where $$Z_{N,\beta} = \sum_\sigma e^{-\beta H(\sigma)}$$ is the partition function, and Onsager's explicit expression for the thermodynamic limit of $$\frac{1}{N}\log Z_{N,\beta}$$ is given in the links above.
• I am using Python. It used to be faster but I started using the decimal library as the partition function grows too large for larger systems for the inbuilt float type. Calculations with the decimal type seem to really slow down the program. As for your suggestion about specific heat, I don't really see how estimating the specific heat would be easier when it is equal up to a scalar constant to $C_{v}$. Unless you are suggesting that I should use the exact expression? That would solve my current problem but I will need to move on to more complicated systems later so it isn't really an option. Commented Aug 24, 2020 at 14:40
• By measure the variance directly, do you mean use the formula $Var(E)=\sum_{i}{p(E_{i})(E_{i}-\left<E\right>)^{2}}$? Not really sure how that would help though? I'll try it out though. Thanks. Commented Aug 24, 2020 at 15:50
• I mean computing the sample variance. The above should help, because your formulas above make ugly factors of $T$ (or $k_B$) appear, which seemed to be your problem, right? In what I suggest, the only relevant thermodynamic parameter is $\beta$ (as it should). Commented Aug 24, 2020 at 16:00
• Yes, as this only differs from the usual definition by a multiplicative constant. (Or simply use $\sigma_E^2$ as it contains the same information, and also allows you a comparison with theoretical results. But I am a mathematical physicist and so couldn't care less about all these constants: I never distinguish between $\beta$ and $1/T$ in the first place. ;) ) Commented Aug 25, 2020 at 5:40