Mean Field - Statmech4ChemBio

MFA applied to the Ising Model¶

In the Ising model, the mean-field approximation (MFA) decouples spin–spin interactions, reducing the system to a single macroscopic parameter: the average magnetization per spin, denoted by $m$ .

\langle s_i s_j \rangle \approx \langle s_i \rangle \langle s_j \rangle = m^2

(2)

\langle H \rangle \approx -J \sum_{\langle ij \rangle} \langle s_i \rangle \langle s_j \rangle - h \sum_i \langle s_i \rangle = -\frac{1}{2} J q m^2 - h m

(3)

The factor of $\frac{1}{2}$ avoids double-counting spin pairs, and $q$ is the coordination number—for instance, $q = 4$ for a 2D square lattice and $q = 6$ for a 3D cubic lattice.

Derivation of the Mean-Field Hamiltonian

Consider the 2D Ising model where each spin $s_i$ interacts with its $q$ nearest neighbors:

H_i = -J \sum_{\delta} s_i s_{i+\delta} - h s_i = -\left( J \sum_{\delta} s_{i+\delta} + h \right) s_i

(4)

Define the local mean field at site $i$ as:

H_i^{\text{MF}} = J \sum_{\delta} s_{i+\delta} + h

(5)

Under the mean-field approximation, replace neighboring spins with their average value:

H_i^{\text{MF}} \approx J q m + h

(6)

The total Hamiltonian then becomes:

H \approx - (J q m + h) \sum_i s_i

(7)

Taking the thermal average and accounting for double counting leads to the previously derived mean-field expression for $\langle H \rangle$ .

Self-Consistent Equation for Magnetization

Partition Function: After applying MFA, the system becomes effectively non-interacting, and the total partition function factorizes as $Z = z^N$ , where $z$ is the partition function of a single spin. For $s_i = \pm 1$ :

E(s_i = \pm 1) = \mp (J q m + h)

(8)

z = e^{\beta (J q m + h)} + e^{-\beta (J q m + h)}

(9)

The probability that a spin is in state $s_i$ is:

p(s_i) = \frac{e^{\beta s_i (J q m + h)}}{z}

(10)

Average Magnetization: The mean value of a spin is then:

m = \langle s_i \rangle = \sum_{s_i} s_i \, p(s_i) = \frac{e^{\beta (J q m + h)} - e^{-\beta (J q m + h)}}{e^{\beta (J q m + h)} + e^{-\beta (J q m + h)}}

(11)

m = \tanh(\beta (J q m + h))

(12)

This self-consistent equation relates the average magnetization $m$ to itself via the hyperbolic tangent function, capturing the thermal competition between spin alignment and disorder.

Free Energy of Mean-Field Ising Models¶

Defining the Macrostate $M$

Consider a spin lattice of $N$ sites, where each spin can point up or down. The total magnetization is given by:

N = N_{+} + N_{-}, \qquad M = N_{+} - N_{-}

(13)

Defining the magnetization per spin $m = M/N$ , the probabilities of spin-up and spin-down states are:

p_{\pm} = \frac{1 \pm m}{2}

(14)

Entropy, Energy, and Free Energy

Entropy is given by the Shannon expression using the probabilities of spin states:

S(m) = -N k_B \left[ \left( \frac{1+m}{2} \right) \log \left( \frac{1+m}{2} \right) + \left( \frac{1-m}{2} \right) \log \left( \frac{1-m}{2} \right) \right]

(15)

Energy in the mean-field approximation is the sum over spins in the presence of an average field:

E(m) = -N \left[ \frac{Jq}{2} m^2 + h m \right]

(16)

Free Energy is the usual Legendre transform:

F(m) = E(m) - T S(m)

(17)

Dimensionless Free Energy per Spin (useful for analysis and plotting):

f(m) = \frac{F}{N k_B T} = - \left[ \frac{Jq}{2 k_B T} m^2 + \frac{h}{k_B T} m \right] + \left[ \left( \frac{1+m}{2} \right) \log \left( \frac{1+m}{2} \right) + \left( \frac{1-m}{2} \right) \log \left( \frac{1-m}{2} \right) \right]

(18)

Finding the Critical Temperature $T_c$

The critical point is determined by the second derivative of the free energy:

\frac{\partial^2 f}{\partial m^2} \bigg|_{m=0} = 1 - \frac{Jq}{k_B T} = 0

(19)

\Rightarrow \quad T_c = \frac{Jq}{k_B}

(20)

Magnetization as a Function of Temperature

The equilibrium magnetization minimizes the free energy, leading to a self-consistency condition:

m = \tanh \left( \frac{Jq}{k_B T} m + \frac{h}{k_B T} \right)

(21)

Rearranged as:

- \left( \frac{Jq}{k_B T} \right) m - \frac{h}{k_B T} + \tanh^{-1}(m) = 0

(22)

For zero external field ( $h = 0$ ), the equation simplifies to:

m = \tanh \left( \frac{T_c}{T} m \right)

(23)

For small $m$ , expand $\tanh^{-1}(m) \approx m + \frac{1}{3} m^3$ , yielding:

- \frac{T_c}{T} m + m + \frac{1}{3} m^3 = 0

(24)

Solving this, we find the magnetization near the critical point behaves as:

m = \pm \sqrt{ \frac{3(T_c - T)}{T_c} }

(25)

Visualizing MF predictions¶

The $h=0$ MFA case

The equation can be solved in a self-consistent manner or graphically by finding intersection between:

$m =tanh(x)$
$x = \frac{Jqm}{k_BT}$

When the slope is equal to one it provides a dividing line between two behaviours.

k_B T_c =qJ

(26)

m = tanh \Big(\frac{Tc}{T} m \Big)

(27)

MFA shows phase-transition!

import ipywidgets as widgets
from ipywidgets import interactive, interact
import matplotlib.pyplot as plt
import numpy as np

# Constants
J = 1  # Interaction strength
k = 1  # Boltzmann constant (set to 1 for simplicity)

def entropy(m):
    # Handle log of zero by adding a small number to the argument
    small_number = 1e-10
    return -(1+m)/2 * np.log((1+m)/2 + small_number) - (1-m)/2 * np.log((1-m)/2 + small_number)

def free_energy(m, T):
    # Calculate the free energy for given m and T
    return -0.5 * J * m**2 - T * entropy(m)

def plot_free_energy(T):
    # Magnetization range
    m_values = np.linspace(-1, 1, 100)
    F_values = [free_energy(m, T) for m in m_values]
    
    # Plotting
    plt.figure(figsize=(8, 5))
    plt.plot(m_values, F_values, label=f'T = {T}')
    plt.xlabel('Magnetization (m)')
    plt.ylabel('Free Energy per Spin (F)')
    plt.title('Mean-field Free Energy vs. Magnetization')
    plt.legend()
    plt.grid(True)
    plt.show()

interactive(plot_free_energy, T=(0.1, 2, 0.1 ))

def mfa_ising_Tc(T=1, Tc=1):

    x  = np.linspace(-3,3,1000)
    
    f = lambda x: (T/Tc)*x
    m = lambda x: np.tanh(x)
    
    plt.plot(x,m(x), lw=3, alpha=0.9, color='green')
    
    plt.plot(x,f(x),'--',color='black')
    idx = np.argwhere(np.diff(np.sign(m(x) - f(x))))
    plt.plot(x[idx], f(x)[idx], 'ro')
    
    plt.legend(['m=tanh(x)', 'x'])
    
    plt.ylim(-2,2)
    plt.grid('True')
    plt.xlabel('m',fontsize=16)
    plt.ylabel(r'$tanh (\frac{Tc}{T} m )$')
    plt.show()

import plotly.graph_objects as go
import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import root_scalar  # Importing root_scalar

def compute_xcross(T, h_over_T):
    
    def f(M):
        Tc = 1.0  # Normalized temperature
        h = h_over_T * T  # Compute actual h from h/T and T
        return M - np.tanh((M + h) / (T / Tc))
    
    # Using symmetric interval for root finding to allow negative solutions
    interval = [-2.0, 2.0]

    # Check if a root is likely within the given range
    if f(interval[0]) * f(interval[1]) > 0:
        return 0.0  # If no root likely, return 0
    
    # Find the root using bisection method within the specified range
    result = root_scalar(f, bracket=interval, method='bisect')
    
    return result.root 

# Define ranges for h/T  and T/Tc (from 0.2 to 1.5)
h_over_T_values = np.linspace(-0.2, 0.2, 101)
T_over_Tc_values = np.linspace(0.1, 1.6, 101)

# Create a meshgrid for T/Tc and h/T
H_over_T, T_over_Tc = np.meshgrid(h_over_T_values, T_over_Tc_values)

# Initialize an array for magnetizations
Magnetizations = np.zeros_like(H_over_T)

# Compute M for each (h/T, T/Tc) pair in the meshgrid
for i, T in enumerate(T_over_Tc_values):
    for j, h in enumerate(h_over_T_values):
        
        Magnetizations[i, j] = compute_xcross(T, h)


# Create a 3D surface plot using Plotly
fig = go.Figure(data=[go.Surface(z=Magnetizations, x=H_over_T, y=T_over_Tc, 
                                 colorscale='RdBu', 
                                 contours={
                                     #'z': {'show': True, 'start': -1.0, 'end': 1.0, 'size': 0.1, 'color':'orange'},
                                     'x': {'show': True, 'start': -1, 'end': 1, 'size': 0.05, 'color':'black'},
                                     'y': {'show': True, 'start': 0.5, 'end': 1.5, 'size': 0.1, 'color':'black'}
                                 })])

# Update plot layout
fig.update_layout(
    title='3D Plot of Magnetization M vs normalized temperature and field; T/Tc and h/T',
    scene=dict(
        xaxis_title='h/T',
        yaxis_title='T/Tc',
        zaxis_title='M',
        aspectmode='manual',
        yaxis_autorange='reversed',  # Reverse the Y-axis
        aspectratio=dict(x=2, y=2, z=0.75),  # Adjust these values as needed
        camera=dict(
            eye=dict(x=-4, y=-2, z=2),  # Adjust camera "eye" position for better view
            up=dict(x=0, y=0, z=1)  # Ensures Z is up
        )
    ),
    autosize=False,
    width=900,
    height=900
)

# Display the figure with configuration options
fig.show()

Landau Theory¶

import numpy as np
import matplotlib.pyplot as plt

# ==== Free energy ====
def f(phi, a, b):
    return 0.5 * a * phi**2 + 0.25 * b * phi**4

# ==== Range of phi ====
phi = np.linspace(-2, 2, 500)
b = 1.0  # fixed positive b

# ==== Plot for different 'a' ====
a_values = [1.0, 0.0, -1.0]
labels = ['a > 0 (disordered)', 'a = 0 (critical)', 'a < 0 (ordered)']
colors = ['#1f77b4', '#ff7f0e', '#2ca02c']

plt.figure(figsize=(7, 5))

for a, label, color in zip(a_values, labels, colors):
    plt.plot(phi, f(phi, a, b), label=label, color=color)

plt.axvline(0, color='k', linestyle='--', alpha=0.3)
plt.xlabel(r'$\phi$ (order parameter)')
plt.ylabel(r'$f(\phi)$ (Landau free energy)')
plt.title('Landau Free Energy: Phase Transition Regimes')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# ==== Parameters ====
N = 128             # Grid size
dx = 1.0            # Spatial resolution
dt = 0.01           # Time step
steps = 1000        # Total time steps
save_every = 10     # Frame interval for animation
epsilon = 1.0       # Interfacial width parameter

# ==== Free energy parameters ====
a = -1.0   # linear coefficient (controls stability of φ = 0)
b = 1.0    # non-linear coefficient (typically > 0 to stabilize)

# ==== Landau free energy derivative ====
def f_prime(phi):
    return a * phi + b * phi**3

# ==== Laplacian with periodic boundary ====
def laplacian(field):
    return (
        -4 * field
        + np.roll(field, 1, axis=0) + np.roll(field, -1, axis=0)
        + np.roll(field, 1, axis=1) + np.roll(field, -1, axis=1)
    ) / dx**2

# ==== Initialize field ====
phi = np.random.rand(N, N) * 0.2 - 0.1  # small noise around 0

# ==== Store frames ====
frames = []

# ==== Simulation loop ====
for step in range(steps):
    mu = -epsilon**2 * laplacian(phi) + f_prime(phi)
    phi += dt * laplacian(mu)
    
    if step % save_every == 0:
        frames.append(phi.copy())

# ==== Animation ====
fig, ax = plt.subplots(figsize=(5, 5))
im = ax.imshow(frames[0], cmap='bwr', vmin=-0.05, vmax=0.05)

ax.axis('off')

def update(frame):
    im.set_array(frame)
    return [im]

ani = FuncAnimation(fig, update, frames=frames, interval=50, blit=True)
plt.close()  # avoid double display
HTML(ani.to_jshtml())