8. Interacting systems¶

The previous exercise sets were designed to derive the fundamental equations and concepts in statistical physics. Set 5 focused on the Boltzmann distribution, giving the probability of a microstate, and from the partition function \(Z\) one could derive expectation values of macroscopic observables (such as the average number of bound particles to a substrate, 5.5). These equations allow one to make predictions about large, complex systems, that are in thermodynamic equilibrium, about their mechanical and structural properties, whether they mix or phase separate, and many other macroscopic properties. For example, exercise Section 4.4 revealed why flexible polymers behave like ideal springs. This principle can be recognised in many soft materials around us, including certain types of food, like jellies and pizza dough. Polymers also grant an elastic quality to our own tissue, in a sophisticated way. It is a puzzling thought that the properties of the fantastic variety of materials around us, including organic materials, are somehow encoded in the microscopic properties of the components, and that they can be derived from an appropriate partition function. It is even more puzzling to actually get that information out of the partition function. Formulating the partition function correctly can already be challenging, and solving it is generally very difficult, if not impossible. One may wonder what the use is of a correct equation that contains ‘all the information that one would like to know’ but cannot be solved. It may be similar to having a poem in a foreign script from another culture. That question will be left open for now, although there may be several answers, and we will focus our attention on an interesting system where we can solve the partition function, in an elegant way, using linear algebra. The results are not obvious, and provide information about the behavior of ferromagnets, semi-flexible polymers, random walkers with memory, lipid rafts in cell membranes, and many more seemingly unrelated phenomena.

This set starts with a few exercises about conditional probabilities and correlation functions, and how these functions are essential for obtaining information about structure at the nanoscale, with the application of scattering experiments. The last exercise guides you through the steps, to derive the free energy of the 1-dimensional Ising model, and draw conclusions about its macroscopic properties. It ends with a short Python-code, that you can use to simulate the 2-dimensional Ising model, using a so-called Metropolis Monte-Carlo algorithm.

8.1. Correlations and correlation functions¶

The following relations may be useful:

\[P(A|B)P(B) = P(A, B) = P(B|A)P(A)\]

In this notation \(P(A|B)\) represents a conditional probability, in this case the probability of \(A\), under the condition that \(B\) is given/known. \(P(A,B)\) represents the joint probability of having \(A\) and \(B\) as outcomes. If \(b_1\) and \(b_2\) are the only two possible outcomes of the stochastic variable \(B\), then one can write:

\[P(A) = P(A|B=b_1)P(B=b_1) + P(A|B=b_2)P(B=b_2)\]

A large fraction of the people who make homework tends to pass for a course. For a certain course this fraction is 90%. Of the people who do not make homework, only 50% manages to pass for the exam. Of the entire group, 60% of the people make their homework. A certain student passes for the course, and the lecturer is wondering: “Did this person make homework?” Given only this information, what is that probability?
A biosensor is designed to detect minute concentrations of a certain protein A. If this molecule is captured by the sensor, there is a 99% probability that it will give a signal. Unfortunately, there are also other proteins in the sample that can weakly interact with the sensor. Protein B is also able to cause a signal, but with a probability of 5%. A certain sample contains a concentration of 2 pM of A protein and 1 nM of B protein. When the solution is flushed through the sensor, it is giving a signal. What is the probability that it captured protein A?

8.2. Scattering experiments and correlation functions¶

Expectation values of the form \(\langle \rho (x) \rho(x') \rangle\) can be very meaningful for describing the local structure and dynamics at the nanoscale. Moreover, they provide the link between data from scattering experiments and the actual spatial structure of the material, that scattered the beam (of X-rays, neutrons, electrons,…)

In a uniform, homogeneous material, the density of particles may be the same at any position: \(\langle \rho(x) \rangle \equiv \rho_0\) but the pair correlation \(\langle \rho (x) \rho(x') \rangle\) is not, except for an ideal gas. Every material has a local structure, and the function \(\rho_0 g(\vec{r}) = \langle \rho (\vec{r}) \rho(\vec{r}') \rangle /\rho_0\) can be interpreted as the conditional probability that there is a density \(\rho(\vec{r})\) at \(\vec{r}\), given that there is a particle at \(\vec{r}'\) (and if the average density is the same at every position). So, it describes the local, average density around a particle in the fluid, which may be larger or smaller than the global density. Scattering profiles provide the Fourier-transform of this function. Explaining this would lead to for for this course. It is sufficient to know that gaining insight in the structure at the atomistic scale requires some translation, and that pair-correlation functions provide the key.

Sketch \(g(r)\) for an ideal gas.
Why does \(g(r)\) converge to 1, for large \(r\)?
Sketch \(g(r)\) for sodium ions in water (for concentrations < 0,1 M so that you can use the information of Homework set 7)
Sketch \(g(r)\) for a crystal and liquid water at room temperature. (be welcome to do a quick literature research what this function may look like - there is no way to guess the precise features of this function)

8.3. the Ising model¶

The Ising model was invented by Wilhelm Lenz who gave it to his student Ernst Ising, who solved it in his thesis in 1924. The model was designed to describe ferromagnetism, i.e. the property of certain materials to become magnetised and to retain their magnetisation in the absence of an external magnetic field. Ernst Ising was able to solve the 1-dimensional model, which we are going to do here as well. The 2-dimensional model can also be solved analytically, but is much more difficult (and that is a royal understatement). It was solved by Lars Onsager in 1944 while he was sitting at home, in a time that European universities were closed. Many scholars have held the opinion that he deserved the Nobel prize for it, but he received it for another contribution, the “reciprocal relations” in chemistry. 1 This model, although it was designed to describe something particular in a simplistic way, appeared to contain some very universal features that were recognised in many other systems. What the Ising model means for statistical physics can perhaps be compared to the hydrogen atom in quantum mechanics.

8.3.1. The 1-dimensional Ising model¶

contains \(N\) particles that we will refer to as “spins”, originally representing the magnetic dipole moment of a certain number \(N\) of atoms, pointing in an upward or downward direction. The spins are arranged in a row, and each spin can take two values, so \(\sigma_n = 1\) or \(\sigma_n = -1\), using \(\sigma_n\) as the symbol for the spin of the \(n^\mathrm{th}\) atom. The spin \(\sigma\) does not have a dimension, which makes the algebra a lot cleaner, and also more general. If desired, one could convert the dimension to angular momentum, or magnetic moment. Every spin interacts with an external field \(H\), and with its neighbours by means of a pair interaction \(J\). If the spin is aligned with the external field, and the external field is pointing upward (\(H\) is positive), the energy is \(-H\) and if it is opposite to the external field it is \(H\). In short:

\[E(\{\sigma\}) = - H \sum \limits_{i=1}^N \sigma_i + J \sum \limits_{i=1}^N \sigma_i\sigma_{i+1},\]

where we use that \(\sigma_{N+1}=\sigma_1\) (e.g. the spins are positioned in a circle). This choice makes the algebra easier, and makes every spin equivalent to the others, so that the expectation values of the individual spins are all the same.

What would be a natural choice for a state vector, describing a microstate of the Ising model?
Argue that the energy of a single microstate can be written as:

\[E(\{\sigma\}) = - \frac{1}{2}H \sum \limits_{i=1}^N (\sigma_i+\sigma_{i+1}) + J \sum \limits_{i=1}^N \sigma_i\sigma_{i+1}\]

[it will become clear later, why one would like to write the energy like this…]
What is the probability of a single microstate, given that the system is equilibrated at a constant temperature.
Show that the partition function \(Z\) for this system can be written as:

\[Z = \sum \limits_{\sigma_1 = \pm 1} \sum \limits_{\sigma_2 = \pm 1} ... \sum \limits_{\sigma_N = \pm 1} \prod \limits_{i=1}^N V(\sigma_{i}, \sigma_{i+1})\]
Show that you can write this as:

\[Z = \mathrm{Tr} \left( V^N \right)\]

with \(V(\sigma_i,\sigma_{i+1})\) a function depending on the interaction parameters \(H\) and \(J\).

To reduce the previous expression to this form, it may be useful to interpret \(\sigma_i\) as the row index of the 2 by 2 matrix \(V(\sigma_i, \sigma_{i+1})\).
Find the eigenvalues of \(V\) and using that the trace is basis-independent, show that

\[Z = \mathrm{Tr} \left( V^N \right)= \mathrm{Tr} \left( \Lambda^N \right) = \lambda_1^N + \lambda_2^N \approx \lambda_1^N,\]

where \(\Lambda\) is a diagonal matrix, containing the eigenvalues of \(V\), and \(\lambda_1\) is the largest eigenvalue.
One can now derive an expression for the mean spin. Similar to problem 5.5 one can derive that (no need to do it yourself, but it can be done with the same effort as in 5.5)

\[N \langle \sigma_i \rangle = -\frac{\partial}{\partial H} k_\mathrm{B}T \ln Z\]

Give an expression for \(\langle \sigma_i \rangle\) in terms of the external field \(H\) and interaction parameter \(J\). [The derivation is not too easy. In principle it is only a matter of taking the derivative, but to simplify the final result requires some clever algebraic juggling. Be welcome to just “give” the result. ]

8.4. The spin-spin correlation function in the Ising model¶

With a similar, but more involved derivation as the previous section, one can also derive an expression for the pair-correlation function \(\langle \sigma_i \sigma_{j} \rangle\),

\[\langle \sigma_i \sigma_{j} \rangle - \langle \sigma_i \rangle \langle \sigma_j \rangle = A \left( \frac{\lambda_2}{\lambda_1} \right)^{|j-i|}\]

with the eigenvalues \(\lambda_2 < \lambda_1\) as in the previous section, \(A \equiv \sin^2 2\phi\) and

\[\tan 2\phi \equiv \frac{1}{e^{2\beta J} \sinh (\beta H)}\]

Describe what kind of information the function \(\langle \sigma_i \sigma_{j} \rangle - \langle \sigma_i \rangle \langle \sigma_j \rangle\) provides.
Show that this correlation function drops exponentially as \(e^{- \frac{|i-j|}{\lambda}}\) with a characteristic length scale \(\lambda \equiv [\ln (\lambda_1/\lambda_2) ]^{-1}\).
How does this length scale change with the interaction parameter \(J\) if there is no external field? [so, you may take \(H=0\) for simplicity. Take \(J\) here to be negative.]
What would be the typical structure of your system for large negative values of \(J\)? What for large positive values? (you can conclude this based on conceptual arguments, without doing algebra)
What would a “spin” represent when the Ising model is applied to:
1. a semi-flexible polymer
2. a persistent random walk (a random walk were the probability of the next step depends on what the previous step was)
3. molecules binding cooperatively to a DNA-molecule
For the systems in the previous question, what would \(J\) and \(H\) represent, or be connected to?

8.5. The Ising model in 2 dimensions with a Monte Carlo algorithm¶

For the 1-dimensional Ising model, one can derive an expression for \(Z\), and from that, one can derive useful relations with observables such as the average spin in terms of the interaction parameters, or the ‘domain size’ (how many spins are typically aligned) given by the correlation length \(\lambda\) following from the previous questions.

One can do this as well for the 2-dimensional Ising model, but it requires a truly heroic effort. Few people in the world have gone through this process. For the 3-dimensional Ising model no one has ever found a solution. Part of the reason is that these models can be solved numerically, and that expectation values can be found with stochastic simulation algorithms.

“Monte Carlo algorithm” is a term used for a certain class of simulation algorithms that generate random numbers to evaluate expectation values, or to model stochastic processes, and are named after the city in Monaco famous for its casinos.

8.5.1. Metropolis Monte Carlo¶

Here we will use the “Metropolis Monte Carlo algorithm” (named after a person with the last name Metropolis) to evaluate expectation values for the 2-dimensional Ising model. For this, a Python code is available, where one can adjust the interaction parameters, grid size, temperature, and other physical parameters, to see how the system depends on these parameters. This algorithm generates \(N\) spins in some initial configuration, and then flips one of these spins at random. If this action lowers the energy, the spin is kept in its new configuration. If the energy increases by a value \(\Delta E\), the algorithm generates a random number \(r\) between 0 and 1, based on a uniform probability distribution, and keeps the spin in the new configuration if \(r < e^{-\Delta E}\). Otherwise it will revert the spin to its original configuration. Then, this procedure is repeated a large number \(N\) of times. By doing so, one actually samples the space of microstates according to the Boltzmann distribution, and therefore, this method is able to evaluate expectation values in thermodynamic equilibrium. For the dynamic properties (diffusion, viscosity, …) one would need other methods.

A simple example of a Monte Carlo algorithm is a method that can be used to evaluate the number \(\pi\). This method generates \(N\) sets of two random numbers \(x_i\) and \(y_i\) from a uniform probability distribution between \(-1\) and \(1\), and counts how many times \(N_{\mathrm{in}}\) the numbers obey the following relation \(x_i^2 + y_i^2 \leq 1\). The ratio \(4 N_{\mathrm{in}}/N\) approaches \(\pi\) as \(N\) increases. Explain why this is the case (conceptually). [Hint: making a little drawing in the \(x,y\)-plane may be helpful]
Use the python code to see what the 2-dimensional Ising model typically looks like:
1. at high and low temperature
2. for large and small negative values of \(J\) and \(H\)
3. for positive values of \(J\)
Did the results of the previous question make sense, physically?
Bonus: Find the critical point. How did you find it? [hint: trial and error is perfectly acceptable - but one could use several forms of information to check whether the system is close to the critical point or not]

1: Supposedly, Onsager’s wife asked the Nobel prize committee, after they congratulated her on the phone, whether it was the Nobel prize in Physics or in Chemistry. A similar confusion arose at the moment that Onsager got a position at Yale, when they discovered he did not have a PhD. He could have used one of his papers as a dissertation, but decided to work on a new topic, and wrote a thesis on “Solutions of the Mathieu equation of period \(4\pi\) and certain related functions”. The physics and chemistry department did not know how to decypher it, but then the mathematics department stepped in, and said that if chemistry would not award the PhD, they would. Onsager became professor before his dissertation was finished.

NB2220 Statistical Physics