# Calculating the temperature in Creutz algorithm

Creutz algorithm is, as far as I know, a less popular Monte Carlo algorithm than the Metropolis algorithm - both are statistical mechanics algorithms.

In one of his papers, Creutz shows his algorithm applied on the Ising model. In the Creutz algorithm a "demon" is used, which can take or give energy to the system. The distribution of the energy of this demon follows the Boltzmann distribution.

Creutz argues that it is "easy" to get the temperature: just calculate the average energy of the demon and get the temperature from there. However, I don't see how this is easy...

# The problem: Calculating $\beta$ from the Boltzmann distribution

The average energy is given by:

$$<E>=\sum_s E_sP(E_s)$$

Because this follows the Boltzmann distribution we have:

$$P(E_s) = \frac{1}{Z}\exp(\frac{-E}{kT})$$

Where Z is the normalization factor. Plugging this back in the first equation:

$$<E> = \frac{\sum_s E_s \exp(\frac{-E_s}{kT})}{\sum_s \exp(\frac{-E_s}{kT})}$$

I don't see how this is easy to invert in order to get T. If there are two energies available I cannot find an analytic expression.

# Creutz expression for $\beta$

In the paper describing his algorithm, Creutz derives an expression: (equation 8 in the paper)

$$\beta = \frac{1}{4} \log(1+\frac{4}{<E>})$$

$\beta$ is the inverse temperature. Also, quoting Creutz: "For a continuous system where the energy can take any positive value, this simplifies to:" (equation 9 in the paper)

$$\beta = \frac{1}{<E_d>}$$

Here $E_d$ is the energy of the demon.

I don't understand how this expression is derived. I have the feeling the equipartition theorem is used.

Creutz has also written a paper on the 3D Ising model where he has also included the Fortran source code of the program. The BETA function calculates the inverse temperature from an average energy E. (Note that this function might only be valid for the 3D Ising model)

# Questions

Main question: How do I get $\beta$ from a Boltzmann distribition?

Follow up question about implementation: If I have multiple demons and want to calculate the temperature of the system, should I take the distributions (over "time") of all demons?

• Good grief that code is horrible. Neglecting the fact that it's FORTRAN77 (since it was from '83), you have unoptimized statement functions and bisection method to find the root of a function. – Kyle Kanos Aug 10 '17 at 13:09

$$P(E) = \frac{\exp(-\beta E)}{Z(\beta)},$$ where $Z(\beta)$ is a normalization factor.
If $E$ can be any positive real number, this is the exponential distribution with parameter beta, and actually $$P(E) = \beta \exp(-\beta E).$$ The expectation value for an exponential distribution yields: $$\mathbb{E}\left[P(E)\right]=\frac{1}{\beta}$$
To derive Eq.(8), you have to notice that for the Ising Hamiltonian, the demon's energy can take only values that are positive integers and multiple of 4. In other words $$P(E) = \begin{cases} \frac{1}{Z\left(\beta\right)}\exp\left(-\beta E\right) & \text{if }E=4k\text{ with }k\geq0 \\ 0 & \text{otherwise}\end{cases}$$ Then a computation of the expectation value $\langle E\rangle$ of this distribution yields the relationship (8).