# Implementing a Monte Carlo Simulation for the Gaussian Model

I want to implement a Monte Carlo simulation of the 1D Gaussian Model (the continuous generalisation of the Ising Model). That is the statistical mechanical model with the following Hamiltonian:

$$H = -\frac{1}{2} \sum_{i=1}^N q_i^2 + K \sum_{i=1}^N q_i q_{i+1}$$

where $$q_i \in \mathbb{R}$$. I have not been able to find much information on this model and it may go by another name I am not aware of (I don't see why it is called the Gaussian Model). Since the spins $$q_i$$ are real valued, naively transferring the Ising Monte Carlo for the Ising Model results in the following difficulty: a random state must be selected but since the spins are real valued there is no clear way to do this. Am I correct in assuming the simulation will converge to the correct result if they are states are drawn from any multivariate distribution which may return any value on the entire real line? Specifically a multivariate Gaussian with no correlation between spins, mean zero and standard deviation one, and then carrying out the accept-reject step as usual?

The Gaussian model was introduced "along the way" by TH Berlin and M Kac, Phys Rev, 86, 821 (1952) as a simplified version of the spherical model. They give a detailed description there, including an exact solution for 1D, 2D and 3D cubic lattices. I've also seen the exact solution in The Critical Point: A Historical Introduction To The Modern Theory Of Critical Phenomena by C Domb, CRC Press (1996). I don't think that I've seen much discussion of the Gaussian model in standard texts on statistical mechanics, possibly because it has a problem at low temperatures (see below).

It is called the Gaussian model because the Boltzmann factor takes the form of a multivariate Gaussian in the spin variables $$\mathbf{q}=\{q_1,q_2,\ldots,q_N\}$$, so the partition function may be written $$Q_N \propto \int_{-\infty}^\infty dq_1 \cdots \int_{-\infty}^\infty dq_N \exp\left[ -\frac{1}{2}\sum_j q_j^2 + K \sum_{\langle j,k\rangle}q_jq_k\right] .$$ The second sum $$\sum_{\langle j,k\rangle}$$ is over nearest neighbour pairs, counting each pair once; for a 1D system the quantity in square brackets is identical with your Hamiltonian. Take care, though, that a factor of $$k_BT$$ (where $$k_B$$ is Boltzmann's constant and $$T$$ the temperature) has been included in the definition of $$K=J/k_BT$$, where the interaction energy between neighbouring spins is $$-Jq_jq_k$$. So, as $$T$$ varies, $$K$$ varies, but the preceding term $$-\frac{1}{2}\sum_j q_j^2$$ does not. Also take care that a negative sign (which we usually put in the Boltzmann formula $$\exp[-\mathcal{H}/k_BT]$$) has also been absorbed into the above formula. So, it can be misleading calling it a Hamiltonian! I would write the Hamiltonian corresponding to that partition function as $$\mathcal{H} = +\frac{1}{2}k_BT\sum_j q_j^2 -J \sum_{\langle j,k\rangle}q_jq_k$$ (note the sign of the first term) but let's stick to your nomenclature $$H=-\mathcal{H}/k_BT$$ (which more-or-less matches Berlin and Kac's formulae).

So the exponent is a quadratic form $$\mathbf{q}\cdot\mathbf{A}\cdot\mathbf{q}$$, where $$\mathbf{A}$$ is an $$N\times N$$ matrix, and the integral over the Gaussian can be evaluated exactly, provided we can calculate the eigenvalues of $$\mathbf{A}$$. This is fairly straightforward (Berlin and Kac show how to do it) for a cubic lattice, a square lattice, and, in particular, for a linear chain, if you assume periodic (wraparound) boundaries. I'll just give the result, eqn (8a) of Berlin and Kac, for the free energy per particle in 1D in the thermodynamic limit $$N\rightarrow\infty$$ $$-\frac{F}{k_BT} = \frac{\ln Q_N}{N} = -\frac{1}{2} \frac{1}{2\pi} \int_0^{2\pi} d\omega \ln[1-2K\cos\omega] .$$ [Take care again: Berlin and Kac confusingly have $$4K$$ here, not $$2K$$: they define $$K=J/2k_BT$$, because their sum over pairs counts every pair twice. I don't do that.]

Here we see the deficiency of the model: for $$K>K_c=1/2$$ the argument of the logarithm becomes negative. At temperatures $$k_BT/J the model breaks down. I've seen $$T_c$$ described as a "critical point" (because quantities diverge there) but it is not a normal phase transition, since there is no thermodynamic state defined for $$T. Don't ask me what happens if you try to simulate the system below this temperature! I don't know!

Finally, I'll turn to your question about sampling in Monte Carlo. The usual approach is to choose one spin at random, define a trial move which consists of updating that spin to give a new state of the system, and accept or reject that move on the basis of the energy change associated with that spin change. If you reject, you keep the spin at its original value (and count this state again, in the computation of simulation averages).

The general approach for continuous degrees of freedom is to update a spin according to $$q_j'= q_j+\Delta q$$ where $$\Delta q$$ is a uniformly sampled random number on some predefined range $$[-\Delta q_{\text{max}},+\Delta q_{\text{max}}]$$. The range is chosen empirically, to give a reasonable acceptance rate of trial moves; it is kept fixed throughout the simulation. It would also be fine to sample $$\Delta q$$ from a Gaussian distribution (symmetric about $$0$$) if you prefer, and you can choose the width of that distribution, again, to give a reasonable acceptance rate; again, keep it fixed throughout. The essential aspect is to satisfy detailed balance: the probability density for choosing the forward move must equal that for the (hypothetical) reverse move that would take $$q_j'$$ back to $$q_j$$. Then the acceptance probability will be given by the Metropolis formula based on the change in the total hamiltonian resulting from the change in $$q_j$$: $$\Delta H =\Delta H_1 + \Delta H_2 = \tfrac{1}{2}({q_j'}^2-q_j^2) + K(q_j'-q_j)(q_{j-1}+q_{j+1}).$$

In this particular model, it is possible to sample in a different way which is quite close to what you say in your question. In other words, having chosen a spin $$j$$ at random, choose $$q_j'$$ afresh from the Gaussian distribution corresponding to the single-spin term in the Boltzmann factor, namely $$\exp[ -\frac{1}{2} {q_j'}^2]$$ (suitably normalized). If you do this, then the acceptance criterion should be the Metropolis formula based only on the change in the neighbour interactions, i.e. just $$\Delta H_2 = K(q_j'-q_j)(q_{j-1}+q_{j+1}).$$ Attempting a complete update of all the spins, and accepting or rejecting that move on the basis of the resulting energy change, is possible in principle, but generally not recommended as the acceptance probability tends to be very low (unless the spins hardly change at all).

Once more, a final warning about those signs. We usually write the Metropolis formula as $$\text{Prob(accept)} = \min \left(1,\exp[-\Delta \mathcal{H}/k_BT] \right)$$ but in your case I guess this should be written $$\text{Prob(accept)} = \min \left(1,\exp[\Delta H] \right)$$ I recommend that you check these details yourself!