Averages of absolute values in Monte Carlo simulation of Ising Model

Question

Consider the 2D Ising model in $0$ field, with Hamiltonian

$$ H=J\sum_{\langle i,j\rangle}\sigma_i\sigma_j$$

The magnetization per spin is defined as

$$M=\frac{1}{N}\sum_i \sigma_i $$

Where $N$ is the number of spins. Since $H$ is symmetric under the flipping of every spins, it seems to me that for all $i$, $\langle\sigma_i\rangle=\langle -\sigma_i\rangle=0$. Put differently, for every configuration with magnetization $M$ there is a configuration with magnetization $-M$ with equal probability, hence $\langle M \rangle=0$. Wikipedia claims that this is a "bogus" argument that only works on a finite volume, and that is fine by me, but it sure looks true for finite $N$.

If I understand correctly when doing finite size Monte Carlo simulations, this problem is solved by computing $\langle |M|\rangle$ instead of $\langle M\rangle$. But what if I have another observable $\Phi$ and I want to calculate the covariance

$$\langle M-\langle M\rangle,\Phi-\langle\Phi\rangle\rangle $$

Should I compute this instead?

$$\langle M-\langle |M|\rangle,\Phi-\langle|\Phi|\rangle\rangle $$

In what case should I swap an average for the average of the corresponding absolute value?

Answer 1

The argument is not "bogus", it simply need to be qualified. Indeed, it is at the heart of the phenomena of (Landau) phase transitions, which is built on spontaneous symmetry breaking. You have a system that has some symmetry, yet we claim that it undergoes phase transition which breaks this symmetry - what gives?! In formal terms, we have an analytic expression for the partition function $Z = \sum \exp(-\beta E_n)$ yet we are looking for a non-analytical behavior which characterizes a phase transition.

The answer to these problems is indeed that the analytical nature of $Z$ holds only for finite sums. At the thermodynamic limit we can get non analytical behavior. In terms of symmetry, under ergodicity we assume that the system travels among all possible configurations with probability proportional to $\exp(-\beta E_n)$, but if we now set the magnetization going on one direction, it is clear that there is a huge energy barrier to flipping all the spins in the other direction. When taking the ergodic assumptions we said that averaging over systems is like averaging over long times, but it is clear that for the system to flip all the spins, the time scale is exponentially small in the energy barrier (for a low enough temperature). Even for not very large systems that may be much larger than the age of the universe! So we potentially solved both problems - however it is not yet clear how to address it analytically. As straight-forward averaging over quantities using $Z$ for finite system will still assume ergodicity and the symmetry would be manifest.

So there are two ways. In Monte-Carlo simulations, the system picks one magnetization when lowering the temperature. So each simulation by itself would actually give a finite magnetization that will persist over rather long cycles, at low temperatures. Averaging over different simulations might be problematic as different runs will pick different directions, randomly, so averaging over absolute value should solve that. But you can check that your simulation actually works when you lower the temperature and it randomly choose one magnetization and just sticks with it. It's fun!

A more formal way, which is also an analytical approach to phase transition, is to add a small symmetry-breaking term to the Hamiltonian, and then take a limit of large-system-small-term. So we write $$H = J \sum \sigma_i \sigma_j - h \sum \sigma_i$$ and set $h$ to be very small. At the limit $N\to\infty$, $h\to 0$ such that $Nh\to \infty$, the phase transition will be evident. And again, this works even for not very large systems when doing averages using different numerical methods. For $T>T_c$ the effect of $h$ is small enough and the average magnetization is zero, while for $T<T_c$ the small term $h$ helps the system pick the right direction of magnetization, but the absolute value of the magnetization is not meaningfully affected by it. Again - you can try it at home.

Answer 2

I think there are two different questions in this one.

• First of all, from the point of view of statistical averaging (i.e. from the point of view of the partition function), the states with the opposite magnetization are equally probable and therefore the average magnetization is always zero. However, these opposite magnetization states are separated by a very large energy barrier. Thus, once the system is stuck in one of them, it is very unlikely that it will be found in the other. In experiment or a simulation this is due to the finite observation/simulation time, whereas in the theory of phase transitions the same effect is achieved by taking thermodynamic limit or calculating partial statistical sums.

• The main reason for using $|M|$ in Monte Carlo simulations is that the simulation is typically run multiple times, e.g., for different values of temperature. In every particular instance the simulation is likely to get stuck in a minimum with a particular orientation of the magnetization, but this minimum is not always the same, which creates difficulties when, e.g., plotting the magnetization dependence on temperature. Thus is why the absolute value is used, which allows to see the magnetization being smoothly reduced to zero rather than jump between $|M|$ and $-|M|$. If you are coding MC, it takes little effort to keep track of both $M$ and $|M|$ (and also $M^2 = |M|^2 for calculating the susceptibility).

Answer 3

I think there are some things to mention more clearly:

• This may be obvious, but Monte Carlo simulations do not give any time evolution of your system. If you look at the time series representing the states of your Marcov chain during simulation, it does not represent the real time evolution.
• In very special cases the Marcov chain DOES represent the time evolution of your system (or is at least similar to it). This is the case for the SINGLE-FLIP algorithm used in many Monte Carlo simulations of the ising model. The reason is that also the real system kind of acts in the same way: It flips spins one by one.

Now, to answer your question: The Monte Carlo simulation for the ising model gives you an average magnetization of zero. And this is TRUE also in reality.

I know this seems counterintuitive, but let’s think about it: Have you tried to run a Monte Carlo simulation for a system with more than 1000 sites? You will see that the autocorrelation time gets huge. If you go to macroscopic systems the autocorrelation time tends to be infinity. And as we pointed out above: The simulation (in this case) represents the time evolution of the real system. You will never see a macroscopic ferromagnet switching sign, but only because there are so many sites (the autocorrelation time is large). In fact, if you wait infinitely long, the magnet will switch polarity and in the average the magnetization is zero (as the simulation correctly predicts).

A famous example for this is superparamagnetism. There, you have very tiny ferromagnets and therefore short autocorrelation times. Also in the real world you can observe these small ferromagnets switching polarity. In fact, this is a problem for producing magnetic hard discs because it limits the density of magnetic sites on the disk.

Edit: To clarify the „bogus“ argument of taking an absolute value: Yes this is bogus. You can not simply introduce an absolute value to get the result you want to have. For small systems and large systems the simulation gives zero average magnetization (and this is correct in reality). Only for infinitely large systems the magnetization is not zero (a kind of spontaneous symmetry breaking). However, if you want to calculate the magnetization of a macroscopic ferromagnet at a point in time, you have to take time evolution into account. The macroscopic magnet does stay for a very long time in one of the magnetizations and does not switch between then a lot (as mentioned above).