In [1]:
using Plots, LinearAlgebra

Ising Model 101¶

Key Theoretical Ideas:¶

Universality, Criticality: physical properties of critical point are similar across many systems
characterized by O(n, d): 

    n = dim of order parameter field;
    d = spacetime dimension

exact solution via transfer matrix
low and high temperature expansions
duality, Wilson loops and gauge theory

1D Ising Model: solution via Transfer Matrix¶

consult ch. 13 of Stat. Mech. (3rd edition) by Pathria and Beale.

$$ \begin{aligned} Z(T, B) &= \sum_{\{\sigma\}} \, e^{-\beta H(\sigma)} \\ H(\sigma) &= -J \sum_{i} \, \sigma_i \, \sigma_{i+1} - \mu B \sum_i \, (\sigma_i+\sigma_{i+1})/2 \\ \implies Z &= {\rm Tr} \, \left( T^{N} \right) \end{aligned} $$

where

$$ T = \begin{bmatrix} e^{\beta(J+\mu B)} & e^{-\beta J} \\ e^{-\beta J} & e^{\beta(J - \mu B)} & \\ \end{bmatrix}. $$

We can directly calculate the trace of powers of matrix in Julia. The average spin can also be evaluated:

$$ \begin{aligned} Z &= {\rm Tr} \, \left( T^{N} \right) \\ \langle \sigma \rangle &= {\rm Tr} \, \left( T^N \sigma_z \right)/Z. \end{aligned} $$

But there is an easier way: Note that

${\rm Tr} \, \left( T^N \right) = \lambda_1^N + \lambda_2^N,$ with $\lambda_{1,2}$ are the two eigenvalues of the $T$, i.e.

$$ \begin{aligned} \lambda_{1,2} &= e^{\beta J} \cosh \beta \mu B \times \left( 1 \pm \sqrt{e^{-4 \beta J}+(1-e^{-4 \beta J}) \, \tanh^2 \beta \mu B} \right). \end{aligned} $$

When $N$ is large, the larger eigenvalue dominates, and we arrive at an analytic expression of the partition function: $Z = \lambda_1^N$

From $Z$ we can derive other theromodynamic variables, e.g. the average spin per site can be computed via

$$ \begin{aligned} \langle \sigma(T, B) \rangle &= \frac{1}{\mu \beta N} \, \frac{\partial}{\partial B} \ln Z \\ &\approx \frac{\tanh \beta \mu B}{\sqrt{e^{-4 \beta J}+(1-e^{-4 \beta J}) \, \tanh^2 \beta \mu B}}. \end{aligned} $$

In [2]:
jcoup = 1.0
mu = 1.0
Bnow = 0.01

function sigma(Tnow)

    beta = 1 / Tnow
    K = beta * jcoup
    H = beta * mu * Bnow

    aa = exp(-4 * K)

    function _f(H)
        lam1(H) = cosh(H) * exp(K) * (1 + sqrt(aa + (1 - aa) * tanh(H)^2))

        return log(lam1(H))
    end

    hh = 1E-5
    ret = imag(_f(H + im * hh)) / hh

    # directly
    retf = tanh(H) / sqrt(aa + (1 - aa) * tanh(H)^2)

    return ret, retf
end

T_arr = range(0.05, 1.4, length=25)
scatter(T_arr, tt -> sigma(tt)[1], label = "<sigma>", xlabel = "T")
plot!(T_arr, tt -> sigma(tt)[2], label = "Ref")
Out[2]:
In [3]:
# exact solution using transfer matrix

# don't go crazy with Ns
Ns = 100

jcoup = 1.0
mu = 1.0

Tnow = 0.3
beta = 1 / Tnow

Bnow = 0.01

K = beta * jcoup
H = beta * mu * Bnow

aa = exp(-4 * K)

Tmat = [exp(K + H) exp(-K); exp(-K) exp(K - H)]
sigmaZ = [1 0; 0 -1]

# direct method
Zpart = tr(Tmat^Ns)
spin = tr(Tmat^Ns * sigmaZ) / Zpart

# via eigenvalues
lam_vec = eigvals(Tmat)
ref = lam_vec[1]^Ns + lam_vec[2]^Ns
ref = log(ref) / Ns

# exact 1D Ising results
lam1 = exp(K) * cosh(H) * (1 + sqrt(aa + (1 - aa) * tanh(H)^2))
m_exact = tanh(H) / sqrt(aa + (1 - aa) * tanh(H)^2)

ret1 = [log(Zpart) / Ns, ref, log(lam1), log(2 * cosh(K))]
ret2 = [spin, m_exact]

[ret1, ret2]
Out[3]:
2-element Vector{Vector{Float64}}:
 [3.366702807629811, 3.366702807629811, 3.3666901509897675, 3.3346051580226734]
 [0.996744360763933, 0.9992722465255379]
In [4]:
function Ising1D(Tnow; Ns = 250, jcoup=1.0, mu=1.0, Bnow=0.01)

    # cold start
    latt = ones(Int64, Ns)
    # hot start
    # latt = rand([-1, 1], Ns)

    # 1D Ising Model, MC implementation

    beta = 1.0 / Tnow


    # ancient relics
	#=
    function Ip(i1)
        # handle periodicity
        ret = (i1 + Ns) % Ns
        # map 0 to Ns
        ret += (ret == 0) * Ns
        return ret
    end
	=#

	# now there is a mod1 function
	Ip = i1 -> mod1(i1, Ns)

    function Esite(i1)
        # neighbors
        neig = latt[Ip(i1 - 1)] + latt[Ip(i1 + 1)]
        return -jcoup * neig - mu * Bnow
    end

	function ham()
		ret = 0.0
		for i1 in 1:Ns
			ret += latt[i1] * Esite(i1)/2
		end
		return ret
	end

    mag() = sum([ss for ss in latt]) / Ns

    function sweep!()

        for i1 = 1:Ns

			# compute the change in energy if flipped at i1
			dE = -2 * latt[i1] * Esite(i1)
			
			# accept if dE < 0 or dE is not so large
			if rand() < exp(-beta * dE)
				latt[i1] *= -1
			end
			
		end
    end

    function update!(; N = 40)
        for i1 = 1:N
            sweep!()
        end
    end

    # starts here

    # stir well
    update!(N = 2500)

    mm = 0.0
    Nm = 150
    for i1 = 1:Nm
        mm += mag() / Nm
        update!(N = 25)
    end

    K = beta * jcoup
    H = beta * mu * Bnow
    aa = exp(-4 * K)

    m_exact = tanh(H) / sqrt(aa + (1 - aa) * tanh(H)^2)

    return mm, m_exact

end
Out[4]:
Ising1D (generic function with 1 method)
In [5]:
T_arr = range(0.1, 1.4, length = 40)

m_arr = similar(T_arr)
mex_arr = similar(T_arr)
for i1 in eachindex(T_arr)
    m_arr[i1], mex_arr[i1] = Ising1D(T_arr[i1])
end
plot(T_arr, mex_arr, label = "<sigma> Ising 1D", title = "1D Ising Model")
scatter!(T_arr, m_arr, label = "MC, Ns = 250")
Out[5]:
In [6]:
function Ising2D(Tnow; Ns = 80, jcoup=1.0, mu=1.0, Bnow=0.01)

    # 2D Ising Model, MC implementation

    beta = 1 / Tnow
    # cold start
    latt = ones(Int, (Ns, Ns))
    # hot start
    # latt = rand([1, -1], (Ns, Ns))

    # handle periodicity
	#=
    function Ip(i1)
        ret = (i1 + Ns) % Ns
        # map 0 to Ns
        ret += (ret == 0) * Ns
        return ret
    end
	=#

	Ip = i1 -> mod1(i1, Ns)

    function ham()
        ret = 0.0
        for i1 = 1:Ns, i2 = 1:Ns
            neig = latt[Ip(i1 + 1), i2] + latt[i1, Ip(i2 + 1)]
            ret += (-jcoup * neig - mu * Bnow) * latt[i1, i2]
        end
        return ret
    end

    function Esite(i1, i2)

        # neighbors
        neig = latt[Ip(i1 + 1), i2] + latt[i1, Ip(i2 + 1)]
        neig += latt[Ip(i1 - 1), i2] + latt[i1, Ip(i2 - 1)]

        return -jcoup * neig - mu * Bnow
    end

    # sum over all spins
    spin() = sum([ss for ss in latt])

    function sweep!()

        for i1 = 1:Ns, i2 = 1:Ns

            Eup = Esite(i1, i2)
            Edown = -Eup

            # if it's now up, flip -> dE = -2*Eup
            dE = -2 * latt[i1, i2] * Esite(i1, i2)
   
			# accept if dE < 0 or dE is not so large
			if rand() < exp(-beta * dE)
				latt[i1, i2] *= -1
			end

			#=
			# alt. formulation: heatbath
			pplus = exp(-beta * Eup)
			pplus /= (pplus + 1 / pplus)

			if rand() < pplus
				latt[i1, i2] = 1
			else
				latt[i1, i2] = -1
			end
			=#

        end
    end

    function update!(; N = 20)
        for i1 = 1:N
            lat = sweep!()
        end
    end

    # stir well
    update!(N = 500)

    s1 = 0.0
    s2 = 0.0
    Eav = 0.0

    # difficult to get fluctuation
    Nm = 400
    for i1 = 1:Nm

        tmp = spin()
        s1 += tmp / Nm
        s2 += tmp^2 / Nm

        Eav += ham() / Ns^2 / Nm

        update!()
    end

    mm = s1 / Ns^2
    sus = (s2 - s1^2) / Ns^2 * beta

    return mm, sus, Eav

end
Out[6]:
Ising2D (generic function with 1 method)

Study of Susceptibility (and other Fluctuations)¶

Define $H = \mu B$ to be an effective external field, it is easy to see the average spin per site can be computed via

$$ \begin{aligned} \langle \sigma \rangle = \frac{1}{\beta N_s^d} \frac{\partial}{\partial H} \, \ln Z \end{aligned} $$

It is an intensive observable. It is also an order parameter. But order parameter is not unique, study of fluctuations of order parameter field is also (more) useful. The susceptibility can be defined (on the lattice) as

$$ \begin{aligned} \chi_{\rm lat} &= \frac{1}{\beta N_s^d} \frac{\partial^2}{\partial H^2} \, \ln Z \\ &= \frac{\partial}{\partial H} \, \langle \sigma \rangle \\ &= \beta N_s^d \, \left( \langle \sigma^2 \rangle - \langle \sigma \rangle^2 \right) \\ &= \frac{\beta}{N_s^d} \, \left( \langle (\sum \sigma_i)^2 \rangle - \langle \sum \sigma_i \rangle^2 \right). \end{aligned} $$

  • Do you understand the equivalence of these expressions?
  • Can you keep track of the factor of $N_s$ and $H$?
  • The last one can be evaluated even in the absence of $H$!

This gives us different ways of studying the fluctuations!

Note that

$\chi \neq \frac{\partial}{\partial T} \langle \sigma \rangle$.

Confusingly, the latter is sometimes called the thermal susceptiblity, we shall avoid this naming! As motivated in the study of Mean Field Theory, susceptiblity is related to an inverse mass square

(see here)

$\chi \propto \frac{1}{m^2}.$

In [7]:
Td = 2 / log(sqrt(2) + 1)
T_arr = range(0.5 * Td, 1.6 * Td, length = 26)

m_arr = similar(T_arr)
dm_arr = similar(T_arr)
sus_arr = similar(T_arr)

for i1 in eachindex(T_arr)
    tmp = Ising2D(T_arr[i1], Bnow = 0.01)
    tmp2 = Ising2D(T_arr[i1], Bnow = 0.02)
    m_arr[i1] = tmp[1]
    dm_arr[i1] = (tmp2[1] - tmp[1]) / 0.01
    sus_arr[i1] = tmp[2]
end

plot1 = scatter(
    T_arr / Td,
    m_arr,
    title = "2D Ising Model",
    xlabel = "T/Td",
    label = "<sigma 2D>",
)
plot2 = plot(
    T_arr / Td,
    dm_arr,
    label = "d<sigma>/dH",
    xlabel = "T/Td",
    color = "gray",
    markers = :true,
    line = (dash = :dash),
)
scatter!(T_arr / Td, sus_arr, label = "sus")

plot(plot1, plot2)
Out[7]:
In [8]:
function lnZ(Tnow)
    # Onsager's solution to the 2D Ising Model
    # lnZ per site; see Pathria (2nd ed.) pp. 383
    jcoup = 1.0
    K = jcoup / Tnow
    kappa = 2 * sinh(2 * K) / cosh(2 * K)^2

    dphi = 0.01
    phi_arr = range(0, pi / 2, step = dphi)
    f(phi) = log(1 + sqrt(1 - kappa^2 * sin(phi)^2))
    ret = 1 / pi * dphi * sum([f(phi + dphi / 2) for phi in phi_arr[1:end-1]])
    return ret + log(sqrt(2) * cosh(2 * K))
end

function Eav(Tnow)
    # E = T^2 d/dE lnZ
    # calculate E per site
    ret = Tnow^2 * imag(lnZ(Tnow + im * 0.001)) / 0.001
    return ret
end

Td = 1 / (0.5 * asinh(1))
T_arr = range(0.1, 2.0 * Td, length = 10)
E_arr = [Eav(Tnow) for Tnow in T_arr]
E_lat = [Ising2D(Tnow, Ns = 150)[3] for Tnow in T_arr]

scatter(T_arr, E_lat, xlabel = "T", ylabel="<E>", label = "MC Ns=150")
plot!(T_arr, E_arr, label = "Eav (Onsager)")
Out[8]:

3D Ising Model¶

For dim = 3 (and above), we need something more robust.

Implementation Notes: (useful tools to learn)¶

moving around lattice: ntuple, eachindex, CartesianIndex CartesianIndices
mod1(i1, Ns) to implement periodic boundary condition, with correct 0 -> Ns
optimize grid setup, update (sweep), and measurement
Metropolis, Heatbath, etc.
In [9]:
function Ising3D(Tnow; Ns = 16)

    # 3D Ising Model, MC implementation; 03.2025
    dim = 3
    jcoup, mu, Bnow = 1.0, 1.0, 0.01
    beta = 1 / Tnow

    latt = ones(Int, ntuple(_ -> Ns, dim))
    sites = CartesianIndices(latt)

    # mod1 for CartesianIndices
    function mod1C(x::CartesianIndex)::CartesianIndex
        return CartesianIndex(mod1.(x.I, Ns))
    end

    # Find neighbors
    function findneighbors(x::CartesianIndex)
        neig = 0
        for d = 1:dim
            dx = CartesianIndex(ntuple(i -> (i == d ? 1 : 0), dim))
            neig += latt[mod1C(x + dx)] + latt[mod1C(x - dx)]
        end
        return neig
    end

    # Hamiltonian
    function ham()
        ret = 0.0
        for site in sites
            neig = findneighbors(site)
            ret += -jcoup * latt[site] * neig / 2 - mu * Bnow * latt[site]
        end
        return ret
    end

    # Energy of a site
    function Esite(site)
        neig = findneighbors(site)
        return -jcoup * neig - mu * Bnow
    end

    # Sweep through the lattice
    function sweep(; metro = true)
        for site in sites
            Eup = Esite(site)

            # figure this out: Eup is energy when latt[site] = 1
            dE = -2 * latt[site] * Eup

            if metro
                if dE < 0 || rand() < exp(-beta * dE)
                    latt[site] *= -1
                end
            else
                # Heatbath
                pplus = exp(-beta * Eup)
                pplus /= (pplus + 1 / pplus)
                latt[site] = rand() < pplus ? 1 : -1
            end
        end
    end

    # Update
    function update(; N = 40, metro = true)
        for _ in 1:N
            sweep(metro = metro)
        end
    end

    # Thermalization
    update(N = 150)

    # Measurements
    s1, s2, Eav = 0.0, 0.0, 0.0
    Nm = 100
    for _ in 1:Nm
        tmp = sum(latt)
        s1 += tmp / Nm
        s2 += tmp^2 / Nm
        Eav += ham() / Ns^dim / Nm
        update(N = 10)
    end

    mm = s1 / Ns^dim
    sus = (s2 - s1^2) / Ns^dim * beta

    return mm, sus, Eav
end
Out[9]:
Ising3D (generic function with 1 method)
In [10]:
# value from Creutz et al. (1985) 
# Computer Investigations of the Three-Dimensional Ising Model
# Td = 1/0.2216
Td = 4.513

T_arr = range(0.25 * Td, 1.8 * Td, length = 20)

obs = []

for Tnow in T_arr

    entry = [Tnow, Ising3D(Tnow, Ns = 16)...]
    push!(obs, entry')

end

obs = vcat(obs...)

plot1 = scatter(
    obs[:, 1],
    obs[:, 2],
    title = "3D Ising Model 16^3",
    label = "spin",
    xlabel = "T",
)

plot2 = scatter(obs[:, 1], obs[:, 3], label = "susceptibility", xlabel = "T")

plot3 = scatter(obs[:, 1], obs[:, 4], label = "energy per site", xlabel = "T")

plot(plot1, plot2, plot3)
Out[10]:
In [ ]: