Phase separation on a lattice: Monte Carlo simulations and Flory-Huggins theory

Each lattice site is either solvent or a polymer bead, and each polymer is a self-avoiding chain whose bonds connect orthogonal neighbors. The simulation uses a square lattice with coordination number \(z=4\). Assigning a contact energy \(u_{ij}\) to every nearest-neighbor pair of species \(i\) and \(j\), the total lattice energy is

\[H = u_{pp}\,n_{pp} + u_{ss}\,n_{ss} + u_{ps}\,n_{ps},\]

where \(n_{pp}\), \(n_{ss}\) and \(n_{ps}\) count polymer-polymer, solvent-solvent and polymer-solvent contacts. We adopt the symmetric convention \(u_{pp}=u_{ss}=-\varepsilon\) and \(u_{ps}=0\), so polymers and solvents both prefer like neighbors with the same strength \(\varepsilon>0\):

\[H = -\varepsilon\,(n_{pp}+n_{ss}).\]

Because each polymer (solvent) site has \(z\) neighbors, the lattice constraints \(2 n_{pp} + n_{ps} = z N_p\) and \(2 n_{ss} + n_{ps} = z N_s\) let us eliminate \(n_{ps}\) and \(n_{ss}\). Dropping the configuration-independent terms gives a Hamiltonian that depends only on \(n_{pp}\) and on the combination

\[\chi = \frac{z}{k_B T}\!\left[u_{ps} - \tfrac{1}{2}(u_{pp}+u_{ss})\right] = \frac{z\,\varepsilon}{k_B T},\]

the Flory interaction parameter, which on this \(z=4\) lattice is \(\chi = 4\beta\varepsilon\).

A proposed move is accepted with the Metropolis probability \(\min\{1,\,\exp(-\Delta H/k_B T)\}\). For any single-bead move (reptation, end flip, corner flip, or single-bead displacement), each polymer-polymer contact gained at the new site is matched by exactly one solvent-solvent contact gained at the old site, so

\[\Delta n_{ss} = \Delta n_{pp}, \qquad -\frac{\Delta H}{k_B T} = 2\beta\varepsilon\,\Delta n_{pp} = \tfrac{\chi}{2}\,\Delta n_{pp}.\]

The acceptance probability is therefore

\[\min\!\big\{1,\;\exp\!\big(\tfrac{\chi}{2}\,\Delta n_{pp}\big)\big\},\]

where \(\Delta n_{pp}\) is the change in the number of polymer-polymer contacts. Larger \(\chi\) makes like-like contacts more favorable and drives demixing.

The companion mean-field Flory–Huggins theory uses the same interaction parameter \(\chi\) and writes the free energy per lattice site as

\[g(\phi) \equiv \frac{f}{k_B T} = \frac{\phi}{N}\ln\phi + (1-\phi)\ln(1-\phi) + \chi\,\phi(1-\phi).\]

Its first derivative is

\[g'(\phi)=\frac{\ln\phi+1}{N}-\ln(1-\phi)-1+\chi(1-2\phi),\]

and its second derivative is

\[g''(\phi)=\frac{1}{N\phi}+\frac{1}{1-\phi}-2\chi.\]

The binodal is the coexistence curve: two compositions \(\phi_-\) and \(\phi_+\) lie on it when they share a common tangent to \(g(\phi)\), equivalently when they have equal chemical potential and osmotic pressure. Inside the binodal, phase separation lowers the free energy overall.

The spinodal is defined by \(g''(\phi)=0\). Between the two spinodal branches, \(g''(\phi)<0\), so a uniform mixture is linearly unstable and even infinitesimal concentration fluctuations grow spontaneously. Between the binodal and spinodal the system is only metastable: demixing is favorable, but it requires a sufficiently large fluctuation or nucleus.

The critical point is where the binodal and spinodal meet (\(g''(\phi_c)=0\) and \(g'''(\phi_c)=0\)), giving

\[\phi_c = \frac{1}{1+\sqrt{N}}, \qquad \chi_c = \frac{1}{2}\left(1 + \frac{1}{\sqrt{N}}\right)^2.\]

Physically, \(\chi_c\) is the threshold interaction parameter at which coexistence first appears: below it the free energy has a single minimum, while above it two competing compositions emerge.

You might notice that the behavior of our 2D simulations differs considerably from the mean-field Flory-Huggins prediction; for example, the critical point is noticeably larger than the mean-field prediction. Mean-field Flory–Huggins predicts a critical interaction strength \(\chi_c^{\text{MF}}=\tfrac{1}{2}(1+1/\sqrt{N})^2\). Mapping the \(N=1\) case to the Ising model gives an Ising coupling \(J=\varepsilon/2\), and Onsager’s exact solution then yields \(\beta\varepsilon_c = \ln(1+\sqrt{2})\), i.e. \(\chi_c^{\text{exact, 2D}} = 4\ln(1+\sqrt{2}) \approx 3.53\), well above the mean-field value \(\chi_c^{\text{MF}}=2\). This Onsager critical is shown as a red diamond on the phase diagram when \(N=1\).