M13- Monte-Carlo Simulation of a Phase Transition
Jonathan Fernsler, Dept. of Physics, Cal Poly, San Luis Obispo – modified from “Computational Physics” by Nicholas J. Giordano
We’ve learned about integration using the Monte-Carlo technique. Another use of this method is to simulate a thermodynamical system. A macroscopic object has too many interacting particles (~1023!) to possibly predict the state of all particles at any time, but the average behavior of the system as a whole is predictable. The interacting particles obey the laws of statistical mechanics and when averaged over the macroscopic object, we can extract equilibrium thermodynamics. We will use Monte Carlo methods to simulate a phase transition induced by a change in temperature of a ferromagnet.
Ising Model of Magnetism
Magnetism is an inherently quantum mechanical phenomena, and as Niels Bohr proved, ferromagnetism cannot exist in a classical system. His proof was accomplished before quantum mechanics had been invented, so this was a indication that the world didn’t work according to classical physics. A quantum property of magnetism is the electron’s spin, which has an associated magnetic moment. Ferromagnetism arises when a collection of spins all point in the same direction, which produces a total moment that is macroscopic in size. We would like to understand how the interactions between spins produce ferromagnetism. We also know that when heated, this magnetism is lost as the system undergoes a phase transition (this happens at about 1000K in iron). We would like to understand how this occurs and how magnetic properties depend on temperature.
We will use the Ising model to explore a system of interacting spins to investigate magnetism. We will treat our ferromagnet as a collection of atoms arranged on a lattice with each atom possessing spin = ½ (see figure). To treat the full quantum mechanical interactions between atoms would be very difficult (we would need to include angular momentum and other effects), so instead we will use a greatly simplified model, which nonetheless captures most of the interesting physics. Each spin, represented as an arrow, is able to point along either the +z or –z direction; i.e. up or down. No other orientation is permitted, and we can therefore simplify the ith spin as si = 1 (this is a strange quantum mechanical property that when measuring spin, the wavefunction “collapses” so it is either along the axis of measurement or antiparallel). Let’s initialize a matrix called spins and choose its dimensions:
% define the experimental parameters
clear all; close all;
nx = 10;% number of spins in x direction
ny = 10; % number of spins in y direction
% initialize an nx by ny lattice of randomly generated +1
% or -1 spins
spin = InitSpin(nx,ny);
Each of the so-called Ising spins interacts with other spins in the lattice: in a ferromagnet, this favors parallel alignment of pairs of spins. In the real system, these interactions will be largest between nearest neighbor spins and fall off rapidly for spins farther apart. In our simple model we will only consider interactions between nearest neighbors, so the energy of the system is:
where the sum is over all nearest neighbor pairs <i,j> , and J is the “exchange constant”, which is positive for a ferromagnet. Therefore, the red spin in the figure interacts only with the yellow spins in the above figure. But what do we do if a spin is on the edge of the lattice? We could have several boundary conditions, but we will choose “periodic” boundaries, which will effectively make our lattice seem larger than it really is. A spin on the edge of a boundary will interact with its usual nearest neighbors, but will also interact with the spin on the opposite side of the lattice: i.e. the lattice wraps around on itself. This effectively makes the lattice infinite in size, which is advantageous for simulating the thermodynamics of a macroscopic system. But the system still doesn’t quite behave the way we’d like, because the maximum distance two spins can be away from each other is L/sqrt(2) where L is the size of the lattice. You can compute the energy of the entire spin lattice with the provided function EnergyMC, where we take J=1, which uses functions to take into account the periodic boundary conditions:
energy = EnergyMC(spin);
In our model, two parallel-aligned spins have an interaction energy –J and two antiparallel spins have energy +J. To lower energy the spins will favor parallel alignment. You might think that this would always lead to a set of spins that are always aligned parallel over the whole sample to lower energy. However, we have to consider the disordering effects of temperature. Assuming that our spin system is in equilibrium with a heat bath at temperature T, the probability of find the system in any particular state, , is proportional to the Boltzmann factor
where E is the energy of the state found from our equation for the energy of the system, kB is Boltzmann’s constant, and P is the probability of finding the system in that state. Therefore, low energies E have a high probability of occurring (favoring parallel spins) and higher energies E have a lower, but still finite probability. Each of these states is a particular configuration of spins, which we call a microstate. With a lattice containing N spins, there are 2N possible microstates. We are interested in systems where N is very large (remember N~1023 for a real macroscopic sample!) so the total number of microstates is extremely large. In order to discover the behavior of the system, we need to calculate the probabilities of being these microstates.
We can also measure the magnetization of the system to observe how this changes with temperature. The magnetization is
where is the magnetization of one particular microstate . Now we want to try to calculate the properties of the system and see how they depend on temperature.
Monte Carlo Method
We wish to see how our system interacts with its environment: the heat bath at temperature T. According to statistical mechanics, the role of a heat bath is to exchange energy with the spin system and thereby bring it into equilibrium at some temperature T. As the system gains or loses energy from the bath, spins are flipped causing the system to move to new microstates. Each microstate corresponds to a particular arrangement of spins. The measured value of magnetization depends on probability of find the system in various microstates.
In the Monte Carlo method, you begin with the system in a particular microstate. You then choose a spin (either sequentially or at random) and the energy required to make it flip, Eflip, is calculated. If Eflip is negative (so the energy is lowered by making the spin change direction) then the spin is flipped and the system is in a new microstate. If Eflip is positive (so the energy is raised by flipping the spin) we make a decision. A randomly generated number between 0 and 1 is compared to the Boltzmann factor exp(-Eflip/kBT). If this number is larger than the Boltzman factor, the spin is flipped, if it is smaller the spin stays the same. Hence the system may move to a different microstate, depending on the value of the random number. This algorithm of choosing a spin and deciding whether to flip it is repeated over the whole lattice, which has now entered a new microstate. The function spinMC accomplishes this:
T=1.3;% Define your temperature
spin = spinMC(spin,T); % one cycle of flipping spins
Make sure to look through each of the provided functions to see that you understand how they work and what they are doing.
Now we can explore this system and examine the thermodynamics. Above, we’ve taken the lattice into a new microstate and we can compute the energy of that microstate. However to explore the system more accurately, we want to examine many microstates of the same system to look at averages and fluctuations. This simple model gives a surprisingly accurate picture of magnetization as we’ll see.
You’ll be using lots of for loops to cycle through the lattice and its microstates, so make sure you’re comfortable with that from the previous assignment. Below you will write some functions and a script to explore the system and place your generated plots and observations in a document about the Ising model.
Assignment: M13_MonteCarloPhaseTrans
- Write a function MagMC(spin) which computes the magnetization of the whole lattice. We can use this to watch how the Ising model can melt from a ferromagnetic phase when temperature is increased.
- You now have functions to compute the energy and magnetization and to advance the Ising model to a new microstate. In your script, compute the energy and magnetization per spin for 1000 time steps (cycling through the lattice to find a new microstate) for T = 1.1 and plot them vs. the Monte Carlo time steps. Do these plots make sense? Now compute the average energy <E> and magnetization <Mper spin over all the time steps from these values in your script. Monte Carlo time steps are not actually related to the real time.
- How does the size of the lattice affect your calculations of average energy? Try a range of LxL lattice sizes where L = 2, 4, 6, …20 and find the average energy <E> per spin for 1000 time cycles. Remember that the “correct” thermodynamics is simulated by an infinitely-sized system. This may take a few minutes to find the average magnetization of a very large system. Generate a plot showing how the lattice size affects the calculated energy and explain what is happening.
- We will use a lattice of size 10x10 for the rest of the simulation. Now we want to explore how temperature affects the energy and magnetization. Set up a loop to compute the average energy <E> and magnetization <M> per spin over 1000 time steps and plot these values on two separate plots versus the temperature ranging from T = 0.5 to 2.5 in steps of 0.05. I suggest that you do not reset the system when you go to the next temperature, instead start where you left off (less time steps to equilibrate). There is a phase transition occurring when the temperature changes: explain what is happening to the spin lattice as a function of temperature. I suggest you plot the absolute value of <M>: why would you want to do this?
- You might notice that the data near the phase transition is noisy. If you run your simulation again, there will be significantly different values in this region from one run to the next. This happens because of large fluctuations in the spin orientation as you reach a critical point at the phase transition itself. The specific heat is related to the energy as C = dE>/dT. At a phase transition in an infinitely large system the specific heat diverges at the critical temperature TC, which means that there are very large fluctuations in the spins. We can compute the size of these fluctuations by computing the average energy per spin over 1000 Monte Carlo time steps as we did before and additionally computing the average energy squared <E2> over the same time steps. The variance is defined as . According to fluctuation-dissipation theorem of statistical mechanics, the variance of the energy is related to the specific heat by . Compute the specific heat from this equation for the range of temperatures from problem 4 and plot vs. the temperature. You can take kB = 1 for this plot, as we only want to see how our specific heat changes when we approach the transition. What temperature does the transition occur at and what is happening here?
- Now let’s add a magnetic field to our Ising model where spin up has energy –B and spin down has energy B (e.g. = 1). The interactions betwee spins are still present, of course. You need to add in the change in energy due to the magnetic field into a modified SpinMC function, call it SpinMCB with an input parameter “B” added. We want to see how the magnetization depends on magnetic field. Write a routine to change the magnetic field from 1 to -1, then back to 1 again in steps of 0.05 on a 10x10 matrix computing <M> for 1000 Monte Carlo times steps at each applied field. Don’t reinitialize spins at each new applied field and make sure that this time, you are plotting <M> and not its absolute value. Try it out at T = 1.5 and 1 (above and below the ferromagnetic transition temperature). Below the transition temperature, you should see a cool effect: hysteresis! This is where the system remembers its history, so it the magnetization looks different for increasing versus decreasing fields.
Extra credit. Try repeating the plots for step 4 above for a 3D square lattice. How does this affect the transition temperature? How about the shape of the specific heat near the transition?
M13_MonteCarloPhaseTrans.doc11/3/18Page 1 of 5