Van der Waals’ square gradient theory

OK, I have encountered this theory again, after many years. The idea is to describe separation between two phases as the minimum of a free energy with respect to an order parameter \phi:

g(\phi) = \frac a2 \phi^2 + \frac b4 \phi^2

When becomes negative, the minimum of g changes from 0 to other values, \pm \phi_0. The function then has the celebrated double minimum feature, which features prominently in many symmetry-breaking theories, including of course the appearance of particle mass mass and, you know, the Big Bang and the Universe.

But here we are just considering phase separation in materials. The interface between two coexisting phases must have some associated cost, and the simplest (lowest order) way to include it is by introducing a total free energy functional

F = \int dr \,  f(\phi, \nabla\phi)

f(\phi) = g(\phi) + \frac c2 (\nabla\phi)^2 .

This is also called a London-Linzburg-Landau free energy, also appearing in their theory of superconductors.

Now, parameters ab, and c are not easy to measure (or, at least, estimate) experimentally, but they are related to: the surface tension, the width of the interface, and the magnitude of the bulk equilibrium order parameter (i.e. \phi_0 ). Here I show how to obtain it in a slightly more general setting, since I was not able to find it on the internet (it can be found e.g, in the book by Rowlinson & Widom).

General square-gradient theory

Let us consider a general square-gradient expression

F = \int dr \, f(\phi) \qquad f(\phi) := g(\phi) + \frac c2 (\nabla\phi)^2,

with a g(\phi) which does not need be the previous one.

The usual Euler equations to find an extremum of the free energy functional are

\frac{\delta F}{\delta \phi} = \frac{\partial F}{\partial \phi} - \nabla \cdot \frac{\partial F}{\partial (\nabla \phi)} = 0

This translates into a modified diffusion equation:

\frac{d g}{d \phi } - c \nabla^2 \phi = 0 .

An alternative form of the Euler equations, since the integrand does not depend on space explicitely, is given by the Beltrami identity:

f  - \nabla\phi \cdot \frac{\partial F}{\partial (\nabla \phi)} = 0

This leads us to

g - \frac c2 (\nabla \phi)^2 = G ,

where G is a constant that must be determined. (In fact, this is an alternative Euler expression that applies since the integration variable (the space) does not appear in the integrand of the functional).

Now, let’s consider variations in the x direction only. At the far left the order parameter has value -\phi_0, and at the right, \phi_0. It follows that the space derivative of the order parameter must be zero at these two extremes. This identifies the constant G as g(\phi_0). In other words,

\Delta g -  \frac c2 (\nabla \phi)^2 = 0 ,

where \Delta g = g(\phi) - g(\phi_0) , the excess free energy (we keep using the nabla symbol, but of course it just means a derivative w.r.t. x ).

If we define an excess free energy functional:

\Delta F = \int dr \, \left[  \Delta g(\phi) + \frac c2 (\nabla\phi)^2 \right],

for the equilibrium \phi_\mathrm{eq}(x) profile,

F_0 = A \int dx \, \left[ \Delta g(\phi_\mathrm{eq}) + \frac c2 (\nabla\phi_\mathrm{eq})^2 \right].

By definition the excess free energy of an interface is its area A times its surface tension. Therefore:

\sigma = \int dx \,  \Delta g(\phi_\mathrm{eq}) + \frac c2 (\nabla\phi_\mathrm{eq})^2 ,

when the equilibrium profile is plugged in it.

Again, instead of solving these head on, our previous result yields

\mathrm{(1)}\qquad  \frac c2 (\nabla \phi)^2 = \Delta g   ,

so that the two terms in \Delta F are exactly equal! This permits writing

\sigma = \int dx \,  \left[ c (\nabla\phi)^2 \right],

or also

\sigma = c \int dx \, (\nabla\phi) (\nabla\phi) = c \int dx \,  \frac{d \phi}{dx}  \sqrt{\frac{2 \Delta g}{c}}.

Now, the latter integral really means a change of variables! We may therefore write

$latex  \mathrm{(2)} \qquad  \sigma = \sqrt{ 2 c}  \int_{-\phi_0}^{\phi_0} d\phi \sqrt{\Delta g}.$

This is a very remarkable expression that estates that the surface tension is the area below the square root of the excess free energy function between the two minima. See the Figure for a plot for the LGL, and an interesting numerical value which will serve us later on.


Integral of the square of the double well between its two minima. It seems dumb to apply the square root of a square, but this is intentional, see below. Credits: Wolfram alpha

Notice this form of the surface tension completely circumvents the expression of the profile. A way to obtain it, alternative to solving the diffusion equation, is to use (1)  again, to write

\frac{d\phi}{dx}= \sqrt{\frac{2\Delta g }{c}},

d\phi \sqrt{\frac{c}{2 \Delta g }} = dx,

\mathrm{(3)}\qquad \int_{\phi_1}^\phi d\phi' \sqrt{\frac{c}{2 \Delta g(\phi')}} = x - x_1 .

In the latter, the value of the order parameter must be known at position x_1. I.e. \phi(x_1)=\phi_1. This permits the calculation of the profile by inversion of the resulting x = \phi .

Equations (2) and (3) may be applied to any square-gradient expression, not just the simple LGL simple double well (for example, it can be applied to van der Waals’ most famous expression for liquid-vapour equilibrium).

Application to LGL

Here, we compute this expressions for the simple double-well potential. Let us write again:

g(\phi) =  -\frac{|a|}2 \phi^2 + \frac b4 \phi^4 .

We will only consider the case in which there is phase separation, and a is negative. In what follows, we will just write a for its absolute value.

Now, we define a normalized order parameter \psi=B\phi such that:

g(\phi) = A \left( -\frac 12 \psi^2 + \frac 14 \psi^4 \right) =: A \bar{g}(\psi) .

The idea is that A \bar{g}(x)= - x^2/2 + x^4/4 contains no physical parameters, which are all absorbed in A and B. Equating the two equations,

-\frac a2 \phi^2 + \frac b4 \phi^4  = A \left( - \frac 12 B^2\psi^2 + \frac 14 B^4 \psi^4 \right)  ,

we find



The usefulness of this transformation is more apparent when we use them in the expression for the surface tension. Indeed

\qquad  \sigma = \sqrt{ 2 c}  \int_{-\phi_0}^{\phi_0} d\phi \sqrt{\Delta g} = \sqrt{ 2 c A} \frac1B  \int_{-\psi_0}^{\psi_0} d\phi \sqrt{\Delta \bar{g}} .  

The $latex \sqrt{A}$ appears from the overall prefactor in the energy, while 1/B comes from the change of integration variable.

The last integral contains no parameters whatsoever! We may predict now

\sigma = n \sqrt{ac}\frac ab ,

an expression perhaps more complicated that may have been expected. In it, n is some dimensionless number, very likely not too large or small. By the way, since a is supposed to be proportional to |T-T_0| close to the critical point, this predicts a classical critical exponent of 3/2 for the surface tension. In fact, the minimum stands at \pm \psi_0 = 1. Recalling B, the extremelly famous critical exponent of 1/2 is predicted for the order parameter.

The excess is given by  \Delta \bar{g}= - x^2 /2 + x^4/4 - ( - 1/2 + 1/4). The latter may be written as

\Delta \bar{g}= (1/4) (1- x^2)^2,

an expression in which the double-well feature is quite prominent.

The integral of $latex (1- x^2)$ is computed in the figure: 4/3 (not a hard one since the square root cancels the power of two!). Finally,

\sigma= \frac 13 \sqrt{2ac}\frac ab .

Now, for the profiles. If we include the square gradient term we may define

f(\phi) =  A \left[ \bar{g}(\psi) + \frac 14 L^2 (\nabla\phi)^2 \right] .

The idea here is to capture the typical length scale L of the interface, since the spatial derivative may then be cast as d/d(x/L) (a factor of 2 is introduced in the definition purely for convenience). Therefore:

ABL^2= 2 c \qquad L=\sqrt{\frac{2c}a},

which does not feature b, and predicts a diverging interfacial spacing at the critical point, with a critical exponent of -1/2.

Going back to Eq. (3) we have

\sqrt{\frac{c}{2 A}} \frac 1B \int_{\psi_1}^\psi d\psi' \frac 1{\sqrt{\Delta \bar{g}(\psi')}} = x - x_1 ,

again with A appearing because of the global prefactor, and B from the change of integration variable. In terms of L:

\frac L2 \int_{\psi_1}^\psi d\psi' \frac 1{\sqrt{\Delta \bar{g}(\psi')}} = x - x_1 ,

which makes clear how the length scale is given by L. Now, let us define the van der Waals dividing surface as the point at which the order parameter takes the value of zero, and let us place that surface at the origin. Then,

\frac 12 \int_{0}^\psi d\psi' \frac 1{\sqrt{\Delta \bar{g}(\psi')}} = \frac{x}{L} .


\int_{0}^x d x' \frac 1{1 - x^2 } = \frac 12\log\frac{1+x}{1-x} .

This function is precisely the inverse of the hyperbolic tangent! Therefore we may invert to get

\psi = \tanh\left[\frac xL \right] \qquad \phi = \sqrt{\frac ab} \tanh\left[\frac xL \right] .

Finding out the coefficients

OK, imagine we are given the value of the surface tension, the bulk concentration and the interfacial length. We may write the surface tension as

\sigma = \frac 13 A L,

with the energy density A=a^2/b = a \phi_0^2. Therefore,


From the value of the bulk concentration, we find

b=\frac{a}{\phi_0^2} = \frac{3\sigma}{L\phi_0^4}.

Finally, from the interfacial length we find for c

c= \frac 12 a L^2=\frac 32 \frac{\sigma L}{\phi_0^2}.

For example, Camley and Brown J. Chem. Phys. 135, 225106 (2011), in a study of 2D hydrodynamics, use \sigma= 0.1 pN (units of force because this is actually a line tension in 2D), and an iterfacial width of L=40 nm.

With these numbers, and an order parameter with a value of 1, we would have

a=b=\frac 34 \times 10^{-5} \,\mathrm{J/m}^2 \qquad c=6 \times 10^{-21} \, \mathrm{J}.

A bit crazy on these units, but in more microscopic ones they seem more sensible:

a=b=\frac 3{400} \,\mathrm{pN/nm} \qquad c=6 \, \mathrm{pN}\,\mathrm{nm}.



Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s