Absence of Discontinuity in Specific Heat Plot Simulated by Ising Model

Question

I am working on 2D Ising model, when I plot "specific heat vs temperature", I can't see any discontinuity at critical temperature around Tc ~ 2.7K. I am enclosing results of all other thermodynamic quantities, which are correct except the specific heat.

Edit 1:

As asked by @Daniel Shapero, I ran the simulation for 64 x 64 lattice, results are posted in comments. As you can see the absence discontinuity in both specific heat and susceptibility plots.

Output:

code:

#!/usr/bin/env python
# coding: utf-8

#--------------------------------------------------
## required packages ##
#--------------------------------------------------

import numpy as np
from numpy.random import random
import matplotlib.pyplot as plt

#--------------------------------------------------
## definitions ##
#--------------------------------------------------

def spin_field(rows, cols):
    ''' generates a configuration with spins -1 and +1'''
    return np.random.choice([-1, 1], size = (rows, cols))

def neighbours(x, y, lattice, dim):
    ''' finds the neigbours of a particular lattice point with periodic boundery conditons '''
    top   = (x - 1, y)
    bottom  = ((x + 1) % dim, y)
    left    = (x, y - 1)
    right = (x, (y + 1) % dim)

    return [lattice[left[0], left[1]],
            lattice[right[0], right[1]],
            lattice[top[0], top[1]],
            lattice[bottom[0], bottom[1]]
           ]

def energy_calc(x, y, lattice, dim, J, B):
    ''' calulates the energy of the whole configuration '''
    dE = - J * lattice[x, y] * (np.sum(neighbours(x, y, lattice, dim)) + B)
    return dE

def total_energy(lattice, dim, J, B):
    ''' Energy of the configuration  '''
    TE = 0
    for x in range(dim):
        for y in range(dim):
            TE += - J * lattice[x, y] * (np.sum(neighbours(x, y, lattice, dim)) + B)
    return TE

#--------------------------------------------------
## relavant details about the lattice ##
#--------------------------------------------------

# number of monte carlo simulations
mcs = 100

# sparse averaging
relax_sweeps = 50

# square lattice dimensions
dim = 20

# coupling constant
J = 1.0

# external magnetic field
B = 0.0

# initialisation of lattice with random spins
lattice = spin_field(dim, dim)

n1 = 1.0 / (relax_sweeps * dim * dim)
n2 = 1.0 / (relax_sweeps * relax_sweeps * dim * dim)

#--------------------------------------------------
## main program ##
#--------------------------------------------------

# number of points between temperature range
m = 50

#dummy indices
index = 0

# this array contains thermodynamics quantities
quantities = np.zeros([m,5])

for temp in np.linspace(0.1,5.0,m):
    M1 = M2 = E0 = E1 = E2 = 0
    for sweep in range(0, mcs + relax_sweeps):
        for x in range(0,dim):
            for y in range(0,dim):
                dE = -2 * energy_calc(x, y, lattice, dim, J, B)
                if (dE <= 0):
                    lattice[x,y] *= -1
                elif (np.exp(-1 * dE / temp) >= random()):
                    lattice[x,y] *= -1
                else:
                    continue
        if (sweep > mcs):
            E0 = total_energy(lattice, dim, J, B)
            E1 += E0
            E2 += E0 * E0

    quantities[index,0] = temp
    quantities[index,1] = (n1*E2 - n2*E1*E1) / (temp * temp)
    index += 1

#--------------------------------------------------
## plotting of relevant thermodynamic quantities ##
#--------------------------------------------------

plt.figure(3)
plt.plot(quantities[:,0],quantities[:,1],'k.',label='{0} x {0}'.format(str(dim)))
plt.legend()
plt.xlabel('Temperature')
plt.ylabel('Specific Heat')
plt.title('Ising Model 2D')
plt.show()

```

Answer 1

I sincerely thank @Daniel Shapero for directing me towards this answer.

Discontinuity in the specific heat or susceptibility curves to be visible significantly, you should take much more finer measurements for large number of sweeps, say, I ran the simulation for

1. 1024 steps for system to reach equilibrium
2. 1024 steps for sparse averaging/to measure specific heat finely
3. Temperature ranging from 1.6 - 3.2, with 100 points in-between
4. Lattice I considered is 20 x 20

Output: