Monte Carlo Simulations for 2D Ising Models

Monte Carlo Simulations for phase separation with 2D Ising Models

1. Spin-flip model:

A paramagnetic material undergoes a second-order phase transition to ferromagnetism as the temperature is lowered below the critical temperature. Below the critical temperature, therefore, the spins of all atoms are oriented in the same direction. As a result the material develops non-zero magnetization. This type of transition can be modeled using Monte Carlo. To start with, the lattice in divided into lx*ly squares and each square is assigned a random spin (+1 or -1). The bond energy between the lattice site i and its nearest neighbor j is given by:

(1)

A random atom of the lattice is chosen and the energy change as a result of the reversal of its spin is calculated which is given as:

(2)

where the sum is over all neighbors. If the change in energy is negative, the reversal of spin is accepted. If, on the other hand, the energy change is positive, the spin reversal is accepted with a probability,

(3)

If J is positive (J =1, for our results), then the bond energy for same-spin neighbors is negative whereas the energy for opposite spin neighbors is positive. For an increase in bond energies, the term in the brackets will become less negative for a higher temperature and p will become higher. Thus, more and more flips that increase energy will be selected. Thus, there will be low net magnetization. For a very low temperature, for a spin flip that results in higher energy the term p will be very low and hence, such spin reversals will rarely be selected. Finally, all spins will be either +1 or -1, i.e. all spins will be parallel. Thus, the stable configuration at temperatures below the critical point has all spins parallel. At temperatures above the critical point, however, entropy prevails and a paramagnetic system is observed. The critical temperature for a 2D Ising model is given by:

(4)

In our units, k_B = 1, and J = 1. So, T ~ 2.269. The code for a spin-flip 2D Ising model is given here. The results of the simulation for the flip model with lattice dimension of 100*100 at T = 1.5 are shown in the figures below. The two colors denote different spins.


Figure 1(a). Initial	Figure 1(b). Intermediate	Figure 1(c). 2200 iter

Figure 2 shows the progress of phase separation for T = 1.5 (< Tc = 2.27). Initially the lattice has a random arrangement of spins with equal positive and negative spins. Within very few iterations the spins form 'ferromagnetic domains' and after further iterations, the whole system is at the same spin: either all +1 or all -1. The number of Monte Carlo Steps (MCS) taken for the complete transformation also varies with different runs. For example, the phase separation may occur so that there are two columns (Figure 3), one of +1 spin and the other of -1 spin. In such case, the time taken for all the spins to become parallel may be very large because the driving force is similar for the growth of any of the two domains. This is a problem with the Metropolis algorithm. A remedy may be a block algorithm in which spins of all sites in a particular block are reversed at a time.

Figure 2. 1500 iterations

Figure 3 shows the average magnetization of the 2D lattice as a function of the Monte Carlo Step for several independent runs with lattice dimensions of 100*100 and T = 1.5. As we can see, the time for complete transformation to ferromagnetic (or parallel spin state) varies for different runs. The final magnetization may be positive or negative. Both the cases in which complete transformation took ~100,000 cycles or more involved intermediate states similar to the one shown in Figure 2.

Figure 3. Average magnetization vs MCS for several runs

As expected, simulation at T = 7 showed no phase separation behavior.

2) Vacancy Mediated spinodal decomposition in binary alloys:


Figure 4(a).	Figure 4(b).

Consider a binary mixture of composition A. Random movement of atoms create small fluctuations of composition resulting in two regions with different compositions. For example, a small region with composition A may split into two regions of different compositions (B and C). The total energy of the two phases (regions) lies on the line-segment joining the points B and C. When the curvature of Gibbs energy vs. composition diagram is greater than zero, as in Figure 4 (a), the random fluctuations of composition result in higher energy for the system (Notice that the line segment joining B and C is higher than the Gibbs energy curve). Thus, the stable state of the system is a solution having a uniform composition. On the other hand, when the curvature of the Gibbs free energy vs. composition diagram is negative (Figure 4 (b)), the local fluctuations in concentrations resulting in two regions of different compositions (B and C) cause a decrease in total energy (Notice that the line segment joining B and C is lower than the Gibbs energy curve). Thus the binary mixture is more stable in the two phase form than as a solution of uniform composition A.

When A and B atoms are added to form a mixture of composition some composition x_A, the Gibbs energy of the mixture is given as:

(5)

where the first term on the right is the weighted sum of the initial values, the second term is the change in configuration entropy and the third term is the excess term for regular solutions. A typical phase Gibbs energy diagram for a spinodal decomposition is given in Figure. The negative curvature of the curve in the region bounded by dash-dot lines arises due to the opposite contributions from the excess Gibbs energy of mixing (positive) and the configuration entropy term (negative). However, if the temperature is very high, the absolute value of configurational entropy term is larger than the excess Gibbs energy of mixing for all compositions. Thus, the total Gibbs energy after mixing is less than the total Gibbs energy prior to mixing and the curvature is positive everywhere.

Figure 5. Contributions to Gibbs energy

The code for vacancy assisted spinodal decomposition is given here. A brief description of the algorithm is as follows:

The initial state of the lattice is randomly assigned. A lattice point can have any of the 3 possible values: +1 (A atom) , 0 (Vacant site) and -1 (B atom). Number of vacant sites is very low. It is assumed that the interstitial diffusion is negligible. The lattice is scanned randomly lx*ly times in each Monte Carlo Step. For a vacant site, a random neighbor is chosen and the energy change resulting from exchange with that site is calculated. The exchange is accepted if it results in a lower energy. Otherwise, the exchange is accepted with a probability given by equation 3. Periodic boundary conditions were applied.

Values of the parameters:

J_aa = J_bb = 1;

J_ab = -1;

J_av = J_bv= J_vv= 0;

fraction of vacancies, fV = 0.001

fraction of A atoms, fA = 0.5

The results show typical spinodal behavior at low temperatures for a binary alloy AB with 50-50 composition. Green, white and blue squares correspond to A atoms, B atoms and vacant sites, respectively. The lattice becomes segregated into checkered domains of single species after some iterations. After a very long time the lattice splits into two regions - one with A and the other with B. Note that a periodic boundary condition implies that the atoms at the ends see each other. Therefore, in the final configuration the B atoms have formed a single domain rather than three separate ones.

A simulation for temperature above the critical point showed no phase segregation.

For, T = 1.5, and lattice dimensions of 50*50,


Figure 6(a). Initial	Figure 6(b). 60,000 iter	Figure 6(c). 61,960,000 iter

For T = 7, no phase separation is observed:

Figure 7. 1500000 iter

3. Spin-Exchange model:

This model is similar to the vacancy assisted spinodal decomposition of binary alloys. However, there are only two spins and the exchange of spins between neighbors does not require a third species (like vacancies earlier). A brief description of the algorithm is as follows:

The lattice is divided into lx*ly squares. Values 0 or 1 are assigned randomly to the squares to distinguish between +1 and -1 spins. A random neighbor is chosen for a randomly selected site i. If the two sites have different spins, the energy change caused by the exchange of the spins is calculated. If the energy change is negative, the exchange is accepted, otherwise the exchange is accepted by the probability given by equation 3. The code for this model is given here and the results are shown in the figures below. After some iterations, the model shows phase separation behavior for temperatures below the critical point (Tc ~ 2.27). No phase separation is observed for temperatures above the critical point.

Values of parameters:

J_aa = J_bb = 30;

J_ab = 28;

fraction of +1 spin sites, fp = 0.5

fraction of -1 spin sites, fn = 0.5

For T = 1.5, and lattice dimensions of 100*100,


Figure 8(a). Initial condition	Figure 8(b). 300 iter	Figure 8(c). 129087 iter	Figure 8(d). 3815738 iter

A parameter of interest in such domain growth problems is the pair-wise correlation function. The correlation function Cr(r) was calculated in (1,0) and (0,1) directions as:

(6)

where array b[i][j] is equal to 1 or -1 depending on the nature of the site (i.e. type of atom on the site, or type of spin).

Figure 9. Correlation function vs r for the Vacancy assisted spinodal decomposition model

Figure 10. Correlation function vs r for the spin-exchange model

Figure 9 shows the variation of correlation function with r for different MCS for the vacancy assisted spinodal decomposition model with lattice dimensions of 100*100. Figure 10 is a similar plot for the spin-exchange model the same lattice dimensions. The correlation function is obviously 1 for r = 0. The value of r for which the correlation function becomes negative for the first time is called the correlation length and is a measure of the phase separation. With progress of iterations, the correlation length increases and is bounded by the half-length of the lattice dimension. For an infinite system, the correlation length keeps on increasing. As can be seen from the plot, the increase in correlation length, and hence the phase separation, follows power law.

(7)

Non-linear data fit can be employed to fit the time-correlation length data, and value of the exponent n can be calculated from it. The progress of phase separation is much slower in the spinodal decomposition model because the process is dependent on the vacancy concentration which is quite low.

Link to the codes:

1. Flip.java

2. Vacancy.java

3. Exchang.java

References:

1. S.J. Mitchell and D.P. Landau, Phase Separation in a Compressible 2D Ising Model, Physical Review Letters, 97, 025701 (2006).

2. http://physics.weber.edu/schroeder/software/ (for plotting the lattice in Java).

3. Tao Pang, Introduction to Computational Physics, Cambridge University Press.

4. Mats Hillert, Phase Equilibria, Phase diagrams and Phase Transformations, Cambridge University Press.