Phase transitions are one of the most elegant places where thermodynamics and geometry meet. Landau’s theory gives us a remarkably simple recipe: write down the free energy as a power series in an order parameter, then find the equilibrium state by locating the minima of that free energy. In this article we’ll build a complete, runnable example — the classic second-order (continuous) phase transition — and use Python to both solve it numerically and visualize what’s happening as temperature crosses the critical point.
The Physics: Landau Free Energy
Near a continuous phase transition, Landau proposed expanding the free energy density $F$ as a power series in the order parameter $m$ (for a ferromagnet, $m$ is the magnetization; it could equally be a superconducting gap, a density difference, or any other symmetry-breaking quantity):
$$
F(m, T) = a_0(T - T_c),m^2 + b,m^4 - h,m
$$
where:
- $T_c$ is the critical temperature,
- $a_0 > 0$ and $b > 0$ are material-dependent constants,
- $h$ is an external field conjugate to $m$ (e.g. an applied magnetic field).
Symmetry (there is no reason for $F$ to prefer $+m$ over $-m$ when $h=0$) forbids odd powers like $m$ or $m^3$, which is why only even powers of $m$ appear when $h=0$.
Equilibrium states are the values of $m$ that minimize $F$ at fixed $T$. Setting $\partial F/\partial m = 0$:
$$
2a_0(T - T_c),m + 4b,m^3 - h = 0
$$
For $h = 0$, this factors nicely:
$$
m\Big[2a_0(T-T_c) + 4b,m^2\Big] = 0
$$
which gives two branches:
$$
m_{eq} = 0 \qquad \text{or} \qquad m_{eq} = \pm\sqrt{\dfrac{a_0(T_c - T)}{2b}}
$$
The second branch only exists (is real) when $T < T_c$. This is the mathematical signature of spontaneous symmetry breaking: above $T_c$ the only stable minimum is $m=0$ (disordered phase), while below $T_c$ two symmetric minima appear at $\pm m_{eq}$ (ordered phase), and $m=0$ becomes a local maximum. The order parameter grows continuously from zero as $T$ decreases past $T_c$ — the hallmark of a second-order phase transition.
For a general applied field $h \neq 0$, the cubic equation above has no simple closed form, so we solve it numerically for each temperature — a nice small example of finding free-energy minima computationally rather than analytically.
The Example Problem
We will:
- Fix $a_0$, $b$, and $T_c$, and sweep a small external field $h$.
- For a grid of temperatures $T$, numerically find the free-energy minima $m_{eq}(T)$ by solving $\partial F/\partial m = 0$ and checking $\partial^2 F/\partial m^2 > 0$.
- Plot $F(m)$ at several representative temperatures to see the single-well → double-well transformation.
- Plot the order parameter $m_{eq}(T)$ versus temperature to see the phase transition curve.
- Render a 3D surface of $F(m,T)$ with the minima traced out as a curve on top of it.
Source Code (Google Colab, single cell)
1 | # ============================================================ |
Code Walkthrough
Physical model (free_energy, dF_dm, d2F_dm2). These three functions directly encode the Landau expansion and its first and second derivatives with respect to $m$. Keeping them separate makes the rest of the code read like the math: we minimize $F$ by finding where $dF/dm = 0$, and we confirm it’s a genuine minimum (not a maximum or inflection) by checking $d^2F/dm^2 > 0$.
solve_equilibrium — vectorized minimization instead of a per-point loop. A naive approach would loop over every temperature and call scipy.optimize.minimize_scalar or brentq individually — for a few hundred points this is not slow, but it doesn’t scale well if you want a fine temperature grid or need to sweep multiple parameters. Instead:
- When $h = 0$, the stationary condition factors exactly, so we compute $m_{eq}$ with a single vectorized NumPy expression — no root-finding at all, and it’s exact to machine precision.
- When $h \neq 0$, the stationary condition is a genuine cubic equation. Rather than calling an iterative solver point-by-point, we apply Cardano’s formula directly with NumPy array operations. When the cubic has three real roots (the “double-well” regime), we compute all three candidates at once with
np.stackand pick whichever one actually gives the lowest free energy at that temperature — this is the numerical version of “find the global minimum, not just any stationary point.”
This turns what could be a slow Python loop into O(N) vectorized array math, which is instant even for a temperature grid with tens of thousands of points.
Plot A (fig1) shows $F(m)$ itself at six temperatures. Above $T_c$ you’ll see a single well centered at $m=0$; as $T$ approaches and drops below $T_c$, the curve flattens and then splits into a double well — this is the free-energy landscape literally deforming as the phase transition happens.
Plot B (fig2) is the phase diagram: order parameter versus temperature. Above $T_c$, $m_{eq}=0$ (disordered/paramagnetic phase). Below $T_c$, two branches $\pm m_{eq}(T)$ peel away continuously from zero — the defining shape of a second-order transition (compare this to a first-order transition, where the order parameter would jump discontinuously).
Plot C (fig3) ties both together in 3D: the whole free-energy landscape $F(m,T)$ as a surface, with a red curve tracing the valley floor — literally the path of the equilibrium minima as temperature is dialed down. Watching the single valley bifurcate into two as you scan from high $T$ to low $T$ is the clearest way to see the bifurcation.
Results
=== Landau Free Energy Minimization === a0 = 1.0, b = 1.0, Tc = 1.0, h = 0.0 Order parameter at T=0.0 : m_eq = 0.7071 Order parameter at T=0.5 : m_eq = 0.4994 Order parameter at T=1.5 : m_eq = 0.0000 Done.
Takeaways
The beauty of Landau theory is that it turns a subtle physical phenomenon — spontaneous symmetry breaking at a phase transition — into a purely geometric statement about a function’s minima. Everything about the physics (the existence of a critical temperature, the continuous growth of the order parameter, the symmetric pair of ordered states) falls directly out of watching a single well deform into a double well as one parameter, $T$, is varied. The same computational pattern — write the free energy, differentiate, solve for stationary points, classify them by the second derivative — generalizes directly to first-order transitions (add an $m^3$ term or flip the sign of $b$), tricritical points, and multi-component order parameters, making it one of the most reusable templates in statistical mechanics.
