First order phase transitions and the dynamics of spinodal decomposition

In this note, we address the connection between equilibrium thermodynamics of first order phase transitions and the dynamics of phase decomposition in the spinodal region.

1 Binodal, Spinodal, and the region of phase coexistence

1.1 Basic concepts

A phase transition of first order occurs, when the free-energy of a thermodynamical system F(T,V,N), with temperature T, volume V and particle number N, has two local minima. This means, that the system has two different homogeneous equilibrium states. These minimima are then necessarily divided by an energy barrier, that is, a region of concave F. This has the immediate consequence, that there is a certain region V I < V < V II in which a homogeneous state is energetically disadvantageous compared to a heterogeneous state, where the system’s material is divided into coexisting parts with the densities nI = N∕V I and nII = N∕V II. Realizations of such heterogeneous states are commonly observed, for example, in the liquid/vapour transition of H2O or CO2, where the coexistence of liquid and vapour can readily be prepared experimentally.

In order to illuminate this phenomenon, we will first determine the boundaries of the coexistence region V I,V II and then calculate the free-energy of the heterogeneous state.

In equilibrium, the coexistence of two phases is possible if and only if pressure p, temperature T, and chemical potential μ are equal in both phases, that is, throughout the entire system. Therefore, one can use the condition

μI(T,pI∕II) = μII(T,pI∕II) = const.,

where pI∕II denotes the coexistence pressure, together with the Gibbs-Duhem relation G = jμj, with free enthalpy G, to infer

GI(T,pI∕II) = GII(T,pI∕II)
FI FII = pI∕II(V  − V )
 II    I. (1)
Here, we have used the relation G = F + pV between the thermodynamic potentials. Since we also know that for T = const.
           ∫ VI
FI − FII = −   dV p,

we can combine the equations (1) and (2) to obtain

     dV p = pI∕II (VII − VI).

This result allows to find the coexistence region by a simple geometrical approach: the famous Maxwell construction. Drawing the isotherm of the homogeneous system in a pressure-volume diagram, the coexistence pressure pI∕II corresponds to a horizontal line which is drawn in a way that the areas A and B in Fig. 1 are of equal size. This uniquely defines pI∕II, V I, and V II.


Figure 1: Isothermal lines of CO2 around the liquid-gas phase transition modelled by the van-der-Waals equation of state. The solid lines are, listed from top to bottom, isothermals calculated for T = 320K, T = Tcr = 303.998K, T = 280K, and T = 260K. The binodal is plotted as a dashed line. It delimits the region where homogeneous states are metastable. The spinodal curve is plotted as a dotted line and it marks the boundary of the region where homogeneous states are unstable. Binodal and spinodal meet in at the critical point, which is marked by a black dot.

The energy of the heterogeneous state is given by F = cIFI + cIIFII, where FI = F(T,V I,N) and FII = F(T,V II,N) while cI = NI∕N and cII = NII∕N denote the fractions of molecules in the two phases. Since cI + cII = 1 and V = cIV I + cIIV II, we can deduce the so-called lever rule:

     VII −-V-            -V-−-VI-
cI = VII − VI, and  cII = VII − VI.

Accordingly, we can write the free-energy of the heterogeneous state as

              V   − V      V − V
Fhet(T,V,N ) =-II----FI + ------IFII,
              VII − VI    VII − VI

which is obviously a straight line with the slope ∂Fhet∕∂V = (FII FI)(V II V I), connecting F(T,V I,N) and F(T,V II,N). Since we argued above, that the pressure p = ∂F∕∂V is constant along any isotherm in the coexistence region, we infer that the straight line Fhet is actually tangent to F in both points F(T,V I,N) and F(T,V II,N). It is thus termed the common tangent in the literature. Due to the concavity of the free-energy of the homogeneous state in between the two local minima, the tangent Fhet is lying below F, making the heterogeneous states on the tangent energetically more favourable. This is shown graphically in Fig. 2.


Figure 2: Free energy of CO2 at temperature T = 260K modelled by the van-der-Waals equation of state. In the binodal region, which is shown in light gray, homogeneous states are metastable, since heterogeneous states of lower free energy, lying on the common tangent (dashed line), are possible. The spinodal region is shown in dark gray.

This region in state-space in which the heterogeneous state has a lower overall free-energy is called the binodal region, and it is the region where phase coexistence is possible. Its boundary, the binodal curve, is given by the loci of V I,V II written as functions of another state variable such as the temperature T, yielding curves V I(T),V II(T). If phase coexistence is not possible for all values of T, then the two curves V I(T),V II(T) will meet in a critical point at a temperature Tcr beyond which the phases are no longer distinguishable from each other.

However, the given argument does not claim that a state within the binodal is unstable but only that it is metastable, that is, unstable to finite but not to infinitesimal perturbations. In fact, these metastable states are commonly observed, for example in form of an overheated liquid or an undercooled gas.

So, when does the homogeneous state actually become unstable? Looking at the pressure p = ∂F(T,V,N)∕∂V of the system, we note that any state with ∂p∕∂V = 2F∕∂V 2 > 0, is mechanically unstable, so that the inflection points 2F(T,V,N)∕∂V 2 = 0, when written as functions of another state variable, such as T, demark a boundary, the so-called spinodal curve beyond which no homogeneous state can exist. The region within this boundary, that is the region of negatively curved F(T,V,N), is called the spinodal region. In the next subsection, we will explain this instability by use of an energetic argument.

It is clear that the spinodal region is a subset of the binodal region, since in between every two minima of the free-energy there is necessarily a region of negative curvature. By the same argument, we can deduce that spinodal and binodal meet in the critical point, because the minima smoothly approach each other for T Tcr, so that minima and inflection points finally coincide.

Throughout this thesis, we will describe the surfactant monolayer in terms of its molecular number density n = N∕V . We therefore use the remainder of this subsection to formulate the free-energy using the variable n, that is F˜(T,V,n) := F(T,V,nV ), and write

F˜(T,V,n) = f(T,n)V,

denoting the system’s free-energy density by f(T,n).

The pressure of the system is then given by

p = ∂F-(T,V,-N)-
    ∂V = ∂F˜(T,V,-n)
   ∂V ∂˜F-(T,V,n)
= ∂F˜(T,V, n)
---∂V----- + ∂F˜(T,V,n)
= ∂F˜(T,V, n)
   ∂V + nμ,
where we introduced the chemical potential μ = (1∕V )∂F(T,V,n)∕∂n. Inserting eq. (3), we find
                   ∂f (n )
˜p(T,V,n) = − f(n)+ n-∂n- .

In order to determine the mechanical stability of homogeneous states, we have to calculate

    ∂V = ∂˜p(T,V,n)
   ∂V n-
= ∂2 ˜F(T,V,n)
   ∂V 2 2n-
V(                         )
  1-∂˜F(T,V,n)-− ∂2 ˜F(T,V,n)
  V    ∂V          ∂n∂V
V 2 2 ˜
In case of eq. (3) we obtain
 2 ˜          ˜                2 ˜           2 ˜    2
∂-F2 = 0,   ∂F-=  ∂f(n)V,    -∂-F-=  ∂f,    ∂-F2 = ∂-f2V,
∂V          ∂V     ∂n        ∂n∂V    ∂n     ∂n     ∂n


∂p(T,V,N )  ∂2F (T,V,N )    n2 ∂2 ˜F(T,V,n)
----------= -------2--- = − -2------2----.
   ∂V           ∂V          V     ∂n

We conclude that, no matter if we look at the free-energy in terms of variable volume V or variable density n, the stability is always determined by the second derivative of the corresponding free-energy function.

1.2 An energetic argument on the instability of homogeneous systems within the spinodal region

The conditon of mechanical stability, ∂p(T,V,N)∕∂V < 0, can be substantiated by a simple energetical argument which, despite its illustrativeness, is - at least to the author’s awareness - never given in the literature.

A homogeneous state is unstable when even an infinitesimal perturbation leads to a state of lower free-energy. Let us assume, that an initially homogeneous system of density n0 is perturbed with a perturbation of amplitude dn, so that a volume portion V + of the fixed total volume V 0 now has the density n0 + dn with infinitesimal dn. Material conservation demands, that another portion V 0 V + is perturbed at the same time so that it now has the density n0 V +(V 0 V +)dn. Then the free-energy difference between the perturbed and the initial state is given by

ΔF = V +f(n0 + dn) + (V 0 V +)f( n0V0 − (n0 + dn )V+ )
        0   +V 0f(n0)
= (V + V 0)f(n0) + V +∂f-
|n 0dn + V +1
2 2
+ (V 0 V +)[                                              ]
         ∂f     V+       1 ∂2f   (  V+   )2   2
 f (n0)−  ∂n|n0V-−-V--dn− 2 ∂n2|n0  V--−-V-   dn
               0   +               0    +
+ O(   )
Collecting terms order by order, we see that the dn0 and dn1 contributions exactly cancel each other. The lowest order term of ΔF is therefore given by
       (           )  2           (   )
ΔF =  1  1+ --V+---  ∂-f2|n0dn2 + O dn3 .
      2     V0 − V+  ∂n

Since V + is by definition positive and less than V 0, we conclude that the sign of ΔF for infinitesimal dn is determined solely by the curvature of the free-energy density in the inital state n0. Note, that the above argument holds true also when we assume the perturbation to be limited to a subvolume V ′⊂ V 0. This means, that a local infinitesimal fluctuation around n0 will reduce the free-energy of the system if and only if 2f(n0)∕∂n2 < 0. Thus, we have reproduced the stability conditon given in the previus section by explicitly referring to the change of the system’s free-energy.

2 Dynamics of spinodal decomposition:
The Cahn-Hilliard equation

So far we have only considered equilibrium states and investigated their stability. The actual dynamics of a system within the spinodal region is governed by the Cahn-Hilliard equation [CH58], which describes the time evolution of the spatially varying density ρ(x,t) of the considered thermodynamical system. It can be written as a conservation law

∂tρ(x,t) = − ∇ ⋅J (x,t),

with the flux

J(x,t) = − α ∇δρ(x,t),

where α denotes a mobility factor. In the simple case of constant α, the Cahn-Hilliard equation takes the form

∂tρ(x,t) = α Δδρ(x,t).

It is important to recognize, that the Cahn-Hilliard equation exhibits potential dynamics. This means that an isolated system approaches an equilibrium state, which corresponds to a local minimum of the free-energy functional F[ρ]. In the terminology of nonlinear dynamics this means, that F[ρ] is the Lyapunov functional of the system.

This holds true for arbitrary positive definite mobility factors. To prove this, one has to consider how F[ρ] changes with time. A short calculation yields

tF = V dDx-δF--
δρ(x)tρ(x) ||
|| insert (6)
= V dDx-δF--
δρ(x)∇⋅[        ]
 α∇ -δF--
    δρ(x) |
| int. by parts
= V dDx[  δF     δF ]
 αδρ(x)∇ δρ(x)- V dDxα(   δF  )
  ∇δρ(x)2 ||
|| Gauss’ law
= ∂V d(D1)xα-δF--
δρ(x) V dDxα(       )
  ∇ -δF---
Thus we have decomposed the temporal change of the free-energy into a volume contribution and a boundary contribution. The value of the boundary integral depends, of course, on the specific boundary conditions of the system under consideration. In most relevant cases, such as periodic boundary conditions or also for ρ const. as x →∞, it simply vanishes. For an isolated system, this will always be true, since in general, a system can be regarded as isolated, only if there is no energy flux through the boundaries. This leads us to the result
       ∫       (       )2
∂F  = −   dDx α  ∇ -δF---  ≤ 0,
 t      V          δρ(x)
                  ≥ 0

that is, either ρ is already in an equilibrium state or it is evolving towards one.

Cahn and Hilliard proposed the following free-energy functional for a general isotropic system of nonuniform composition [CH58]:

      ∫     { κ              }
F[ρ] =   dDx  --(∇ ρ)2 + fhom (ρ) .
       V      2

The idea behind this model is simple: To molecules within a subvolume dV of the system it makes little difference whether they are part of a homogeneous system or of an inhomogeneous system with very small density gradients. Thus, as a zeroth approximation, the system’s overall free-energy simply equals the sum of the free-energies of each infinitesimal subvolume, fhom(ρ(⃗x))dV . One can expand the free-energy density in powers of ρ, keeping only the lowest orders:

         2                        (1)        (2)
f(ρ,∇ρ,∇ ρ,...) = fhom(ρ)+ Li∂iρ + κij ∂i∂jρ + κij (∂iρ)(∂jρ)+ ...

The tensors κ(1,2) and the vector L reflect the symmetry properties of the considered material. For isotropic media, one obtains L = 0 and κ(1) = κ1I, κ(2) = κ2I, yielding the free-energy density of eq. (8), where κ = 1∕dρ + κ22.

Inserting the free-energy functional (8) into the Cahn- Hilliard equation (6), we obtain

∂tρ(x,t) = − κΔ2 ρ+ ∂-fhomΔ ρ+ ∂fhom-(∇ ρ)2.
                    ∂ρ2        ∂ρ

From our considerations in section 1 we already know qualitatively what to expect from solutions to this equation: Homogeneous system outside of the binodal region are absolutely stable, within the binodal region they are metastable, and in the spinodal region they are unstable. Upon instability, the homogeneous state decomposes into coexisting domains of different densities ρI, ρII.

Spatially homogeneous fields ρ = ˆρ = const. are always stationary solutions of the Cahn-Hilliard equation. The linear stability of these solutions is readily obtained by looking at small perturbations ζ(x,t). Inserting ρ(x,t) = ˆρ+ ζ(x,t) into eq. (9), and keeping only the linear terms, we find that the time evolution of the perturbation is governed by

       {                }
           2   ∂2fhom-
∂tζ = −  κΔ  −  ∂ρ2  |ρˆΔ   ζ.

Using the ansatz ζ exp(λt + ik x), one finds the dispersion relation of the perturbation, which only depends on the absolute value of the wavevector k = |k|, since only even powers of are applied to ζ:

          {       2      }
λ(k) = − k2 κk2 + ∂-fhom|  .
                   ∂ρ2  ˆρ


Figure 3: Free energy of homogeneous solutions (left) and dispersion relations of small perturbations around these solutions (right). The binodal region is drawn in light gray while the spinodal region is coloured dark gray. For a symmetric Fhom the common tangent (dashed line) is simply the horizontal line connecting the minima of the free energy. Solutions from within the spinodal region, where Fhom is concave, are linearly unstable and have a finite band of unstable modes.

From this expression, it is obvious, that the stability of a homogeneous solution ˆρ is determined by the curvature of the free-energy at ρ = ˆρ . Physically speaking, a solution ρˆ is stable, as long as it is chosen from outside the spinodal region, which is defined as the region of concave fhom. Dispersion relations for three different ˆρ in a system with a free-energy of symmetric double-well shape are shown in Figure 3. Otherwise, there is finite band of modes k with 0 < k < ∘ (∂2fhom-∕∂ρ2)|ˆc which are amplified with λ(k) > 0. The wavenumber kmax, which is maximally amplified, can be calculated from eq. (10) as

       ∘ ------------
           1 ∂2f
kmax =   −------hom2-|ˆρ.
          2κ  ∂ ρ

The corresponding maximal amplification is given by

                  3 ∂2fhom
λmax := λ(kmax) = − 4-∂ρ2-|ˆρ.

With kmax and λmax we have obtained an estimate of the inverse time and length scales of the early stages of spinodal decomposition.

However, after a short time, the nonlinearity of eq. 9 strikes and leads to a different behaviour known as coarsening of the domains which are formed during the initial pattern formation. During this process, large domains grow larger while small domains grow smaller and finally disappear. Two mechanisms can be distinguished: Firstly, two domains can directly merge into a single larger one, thereby minimizing thier overall interfacial energy. Secondly, material diffuses from the surface of small domains towards the larger ones, that is, the domains are absorbed without direct contact to other domains. This phenomenon is known as Ostwald-ripening.


Figure 4: Spacetime plot of spinodal decomposition obtained from direct numerical simulation of the Cahn-Hilliard equation on a periodic domain. Red: ρ = 1, blue: ρ = 1. differences were used for spatial order Parameters: L = 250, N = 512, dt = 4,α = κ = 1.

Figure 4 visualizes the dynamics of this process in a spacetime plot. The domain shape depends on the amount of available material. The phase with the smaller overall volume will form spherical domains in order to diminish the interfacial energy. If, however, both phases have roughly the same total volume, one will rather observe labyrinthine structures. The coarsening process will continue, until finally only one domain of each phase remains, divided by the smallest possible interface.


Figure 5: Solutions of the Cahn-Hilliard equation on a periodic domain, as obtained from direct numerical simulation. Red means ρ = 0, blue means ρ = 1. The initial conditions were randomly perturbed homogeneous solutions. In the first row ˆρ = 0.4, in the second row ˆρ = 0.1, and in the third row ˆρ = 0. In each row, time increases from left to right. The coarsening of the structures towards higher average wavelength can be clearly seen. Comparing the pictures of each column, one can see the different morphologies of the resulting patterns, going from circular domains to labyrinthine patterns as ˆρ 0.


[CH58]   John W. Cahn and John E. Hilliard. Free energy of a nonuniform system. i. interfacial free energy. J. Chem. Phys., 28(2):258–267, 1958.