1. Write a Metropolis Monte Carlo program for simulating a 2-dimensional lattice of N by N Ising spins (you can use any language that you are comfortable with - you just need to ensure that you can create images of the results from your data). For convenience use N=10 (i.e., a 10 by 10 lattice of spins) for initial studies. Use J=1, kB = 1 and h=0 for convenience. For boundaries you can use periodic conditions (i.e., the 2D lattice is actually on a torus.)
There are several nice websites which explain the algorithm (e.g., https://cs.adelaide.edu.au/~paulc/teaching/montecarlo/p_montecarlo.html) but please do not copy any program - you should write it yourself after understanding the algorithm.
It is important to decide what update scheme you will use - try a random updating where each spin whose orientation is to be evolved at a particular time instant is chosen randomly from the N*N spins. If instead you had chosen a parallel update where all spins have their states evolved at the same time you can get stuck in unphysical "checkerboard" patterns (e.g., see the article by Brian Hayes: "The World in a Spin")
2. Starting from randomly oriented Ising spins initially (i.e., with probability 1/2, the initial orientation of a spin is up or down), evolve the system for T=0. Does the system converge to an ordered state at long times ? Show images of the spin orientations on the lattice (e.g., using a suitable color code or arrows as in MATLAB's quiver plot routine) at times t=0, 100, 500 and 5000 Monte Carlo updates (a Monte Carlo update corresponds to N*N successive spin updates).
3. Starting from Ising spins oriented initially in the same direction (all up or all down), evolve the system for T=10. What is the average magnetization per site of the system at long times ? Show images of the spin orientations on the lattice (e.g., using a suitable color code or arrows as in MATLAB's quiver plot routine) at times t=0, 10, 100 and 500 Monte Carlo updates.
4. Calculate the average magnetization per site m of the system from temperatures T=0 to 5, increasing in steps of 0.25. As the initial behavior will depend on the initial state you begin with, you can reject the first 100 MC updates and then average over the next 1000 MC updates.
Plot m as a function of T.
5. Repeat the exercise in (4) but with a larger N=50 lattice. Plot m as a function of T. How does this compare with what you saw for N=10.
Based on this can you venture to say how the m vs T curve will change as N tends to infinity ? Around what value of temperature does the phase transition appear to happen ?