Part III Theoretical Physics of Soft Condensed
Matter
Based on lectures by M. E. Cates
Notes taken by Dexter Chua
Lent 2018
These notes are not endorsed by the lecturers, and I have modified them (often
significantly) after lectures. They are nowhere near accurate representations of what
was actually lectured, and in particular, all errors are almost surely mine.
Soft Condensed Matter refers to liquid crystals, emulsions, molten p olymers and other
microstructured fluids or semi-solid materials. Alongside many high-tech examples,
domestic and biological instances include mayonnaise, toothpaste, engine oil, shaving
cream, and the lubricant that stops our joints scraping together. Their behaviour is
classical (~ = 0) but rarely is it deterministic: thermal noise is generally important.
The basic modelling approach therefore involves continuous classical field theories,
generally with noise so that the equations of motion are stochastic PDEs. The form
of these equations is helpfully constrained by the requirement that the Boltzmann
distribution is regained in the steady state (when this indeed holds, i.e. for systems
in contact with a heat bath but not subject to forcing). Both the dynamical and
steady-state behaviours have a natural expression in terms of path integrals, defined
as weighted sums of trajectories (for dynamics) or configurations (for steady state).
These concepts will be introduced in a relatively informal way, focusing on how they
can b e used for actual calculations.
In many cases mean-field treatments are sufficient, simplifying matters considerably.
But we will also meet examples such as the phase transition from an isotropic fluid
to a ‘smectic liquid crystal’ (a layered state which is periodic, with solid-like order,
in one direction but can flow freely in the other two). Here mean-field theory gets
the wrong answer for the order of the transition, but the right one is found in a
self-consistent treatment that lies one step beyond mean-field (and several steps short
of the renormalization group, whose application to classical field theories is discussed
in other courses but not this one).
Imp ortant models of soft matter include diffusive
φ
4
field theory (‘Model B’), and
the noisy Navier–Stokes equation which describes fluid mechanics at colloidal scales,
where the noise term is responsible for Brownian motion of suspended particles in a
fluid. Coupling these together creates ‘Model H’, a theory that describes the physics of
fluid-fluid mixtures (that is, emulsions). We will explore Model B, and then Model H,
in some depth. We will also explore the continuum theory of nematic liquid crystals,
which spontaneously break rotational but not translational symmetry, focusing on
top ological defects and their associated mathematical structure such as homotopy
classes.
Finally, the course will cover some recent extensions of the same general approach to
systems whose microscopic dynamics does not have time-reversal symmetry, such as
self-propelled colloidal swimmers. These systems do not have a Boltzmann distribution
in steady state; without that constraint, new field theories arise that are the subject of
ongoing research.
Pre-requisites
Knowledge of Statistical Mechanics at an undergraduate level is essential. This course
complements the following Michaelmas Term courses although none are prerequisites:
Statistical Field Theory; Biological Physics and Complex Fluids; Slow Viscous Flow;
Quantum Field Theory.
Contents
0 Introduction
0.1 The physics
0.2 The mathematics
1 Revision of equilibrium statistical physics
1.1 Thermodynamics
1.2 Coarse Graining
2 Mean field theory
2.1 Binary fluids
2.2 Nematic liquid crystals
3 Functional derivatives and integrals
3.1 Functional derivatives
3.2 Functional integrals
4 The variational method
4.1 The variational method
4.2 Smectic liquid crystals
5 Dynamics
5.1 A single particle
5.2 The Fokker–Planck equation
5.3 Field theories
6 Model B
7 Model H
8 Liquid crystals hydrodynamics
8.1 Liquid crystal models
8.2 Coarsening dynamics for nematics
8.3 Topological defects in three dimensions
9 Active Soft Matter
0 Introduction
0.1 The physics
Unsurprisingly, in this course, we are going to study models of soft-condensed
matter. Soft condensed matter of various types are ubiquitous in daily life:
Type Examples
emulsions mayonnaise pharmaceuticals
suspensions toothpaste paints and ceramics
liquid crystals wet soap displays
polymers gum plastics
The key property that makes them “soft” is that they are easy to change in
shape but not volume (except foams). To be precise,
They have a shear modulus
G
of
10
2
–10
7
Pascals (compare with steel,
which has a shear modulus of 10
10
Pascals).
The bulk modulus
K
remains large, with order of magnitude
K
10
10
Pascal. As K/G , this is the same as the object is incompressible.
Soft condensed matter exhibit viscoelasticity, i.e. they have slow response to
a changing condition. Suppose we suddenly apply a force on the material. We
can graph the force and the response in a single graph:
t
σ
0
σ
0
/G
τ
η
1
here the blue, solid line is the force applied and the red, dashed line is the
response. The slope displayed is
η
1
, and
η G
0
τ
is the viscosity. Note that
the time scale for the change is of the order of a few seconds! The reason for
this is large internal length scales.
Thing Length scale
Polymer 100 nm
Colloids 1 µm
Liquid crystal domains 1 µm
These are all much much larger than the length scale of atoms.
Being mathematicians, we want to be able to model such systems. First of
all, observe that quantum fluctuations are negligible. Indeed, the time scale
τ
Q
of quantum fluctuations is given by
~ω
Q
=
~
τ
Q
' k
B
T.
At room temperature, this gives that
τ
Q
10
13
s
, which is much much smaller
than soft matter time scales, which are of the order of seconds and minutes. So
we might as well set ~ = 0.
The course would be short if there were no fluctuations at all. The counterpart
is that thermal fluctuations do matter.
To give an example, suppose we have some hard, spherical colloids suspended
in water, each of radius
a ' 1 µm
. An important quantity that determines the
behaviour of the colloid is the volume fraction
Φ =
4
3
πa
3
N
V
,
where N is the number of colloid particles.
Experimentally, we observe that when Φ
<
0
.
49, then this behaves like fluid,
and the colloids are free to move around.
In this regime, the colloid particles undergo Brownian motion. The time scale of
the motion is determined by the diffusivity constant, which turns out to be
D =
k
B
T
6πη
s
a
,
where
η
s
is the solvent viscosity. Thus, the time
τ
it takes for the particle to
move through a distance of its own radius
a
is given by
a
2
=
Dτ
, which we can
solve to give
τ
a
3
η
s
k
B
T
.
In general, this is much longer than the time scale
τ
Q
of quantum fluctuations,
since a
3
η
S
~.
When Φ > 0.55, then the colloids fall into a crystal structure:
Here the colloids don’t necessarily touch, but there is still resistance to change
in shape due to the entropy changes associated. We can find that the elasticity
is given by
G ' k
B
T
N
V
.
In both cases, we see that the elasticity and time scales are given in terms of
k
B
T
. If we ignore thermal fluctuations, then we have
G
= 0 and
τ
=
, which
is extremely boring, and more importantly, is not how the real world behaves!
0.2 The mathematics
To model the systems, one might begin by looking at the microscopic laws
of physics, and build models out of them. However, this is usually infeasible,
because there are too many atoms and molecules lying around to track them
individually. The solution is to do some coarse-graining of the system. For
example, if we are modelling colloids, we can introduce a function
ψ
(
r
) that
tells us the colloid density near the point
r
. We then look for laws on how this
function
ψ
behave, which is usually phenomenological, i.e. we try to find some
equations that happen to model the real world well, as opposed to deriving these
laws from the underlying microscopic principles. In general,
ψ
will be some sort
of order parameter that describes the substance.
The first thing we want to understand is the equilibrium statistical physics.
This tell us what we expect the field
ψ
to look like after it settles down, and also
crucially, how the field is expected to fluctuate in equilibrium. Ultimately, what
we get out of this is
P
[
ψ
(
r
)], the probability (density) of seeing a particular field
configuration at equilibrium. The simplest way of understanding this is via mean
field theory, which seeks a single field
ψ
that maximizes
P
[
ψ
(
r
)]. However, this
does not take into account fluctuations, A slightly more robust way of dealing
with fluctuations is the variational method, which we will study next.
After understanding the equilibrium statistical physics, we turn to under-
standing the dynamics, namely what happens when we start our system in a
non-equilibrium state. We will be interested in systems that undergo phase
transition. For example, liquid crystals tend to be in a disordered state at high
temperatures and ordered at low temperatures. What we can do is then to
prepare our liquid crystals at high temperature so that it stays at a homogeneous,
disordered state, and then rapidly decrease the temperature. We then expect
the system to evolve towards a ordered state, and we want to understand how it
does so.
The first step to understanding dynamics is to talk about the hydrodynamic
level equations, which are deterministic PDEs for how the system evolves. These
usually look like
˙
ψ(r, t) = ··· ,
These equations come from our understanding of equilibrium statistical mechan-
ics, and in particular that of the free energy functional and chemical potential.
Naturally, the hydrodynamic equations do not take into account fluctuations,
but report the expected evolution in time. These equations are particularly
useful in late time behaviour where the existing movement and changes dominate
over the small fluctuations.
To factor in the fluctuations, we promote our hydrodynamic PDEs to stochas-
tic PDEs of the form
˙
ψ(r, t) = ···+ noise.
Usually, the noise comes from random external influence we do not wish to
explicitly model. For example, suspensions in water are bombarded by the water
molecules all the time, and we model the effect by a noise term. Since the noise is
the contribution of a large number of largely independent factors, it is reasonable
to model it as Gaussian noise.
The mean noise will always be zero, and we must determine the variance. The
key insight is that this random noise is the mechanism by which the Boltzmann
distribution arises. Thus, the probability distribution of the field
ψ
determined
by the stochastic PDE at equilibrium must coincide with what we know from
equilibrium statistical dynamics. Since there is only one parameter we can toggle
for the random noise, this determines it completely. This is the fluctuation
dissipation theorem.
Example.
To model a one-component isothermal fluid such as water, we can
take
ψ
(
r, t
) to consist of the density
ρ
and velocity
v
. The hydrodynamic PDE
is exactly the Navier–Stokes equation. Assuming incompressibility, so that
˙ρ
= 0,
we get
ρ(
˙
v + v · v) = η
2
v p,
We can promote this to a stochastic PDE, which is usually called the Navier–
Stokes–Landau–Lipschitz equation. This is given by
ρ(
˙
v + v · v) = η
2
v p + ·Σ
N
,
The last term is thought of as a noise stress tensor on our fluid, and is conven-
tionally treated as a Gaussian. As mentioned, this is fixed by the fluctuation-
dissipation theorem, and it turns out this is given by
hΣ
N
ij
(r, t
N
k`
(r
0
, t
0
)i = 2k
B
T η(δ
i`
δ
jk
+ δ
ik
δ
j`
)δ(r r
0
)δ(t t
0
).
Example.
If we want to describe a binary fluid, i.e. a mixture of two fluids, we
introduce a further composition function
φ
that describes the (local) proportion
of the fluids present.
If we think about liquid crystals, then we need to add the molecular orienta-
tion.
1 Revision of equilibrium statistical physics
1.1 Thermodynamics
A central concept in statistical physics is entropy.
Definition (Entropy). The entropy of a system is
S = k
B
X
i
p
i
log p
i
,
where
k
B
is Boltzmann’s constant,
i
is a microstate a complete specification
of the microscopics (e.g. the list of all particle coordinates and velocities) and
p
i
is the probability of being in a certain microstate.
The axiom of Gibbs is that a system in thermal equilibrium maximizes
S
subject to applicable constraints.
Example.
In an isolated system, the number of particles
N
, the energy
E
and
the volume
V
are all fixed. Our microstates then range over all microstates
that have this prescribed number of particles, energy and volume only. After
restricting to such states, the only constraint is
X
i
p
i
= 1.
Gibbs says we should maximize
S
. Writing
λ
for the Lagrange multiplier
maintaining this constraint, we require
p
i
S λ
X
i
p
i
!
= 0.
So we find that
k
B
log p
i
+ 1 λ = 0
for all i. Thus, we see that all p
i
are equal.
The above example does not give rise to the Boltzmann distribution, since
our system is completely isolated. In the Boltzmann distribution, instead of
fixing E, we fix the average value of E instead.
Example.
Consider a system of fixed
N, V
in contact with a heat bath. So
E
is no longer fixed, and fluctuates around some average
hEi
=
¯
E
. So we can
apply Gibbs’ principle again, where we now sum over all states of all
E
, with
the restrictions
X
p
i
E
i
=
¯
E,
X
p
i
= 1.
So our equation is
p
i
S λ
I
X
p
i
λ
E
X
p
i
E
i
= 0.
Differentiating this with respect to p
i
, we get
k
B
(log p
i
+ 1) λ
I
λ
E
E
i
= 0.
So it follows that
p
i
=
1
Z
e
βE
i
,
where Z =
P
i
e
βE
i
and β = λ
E
/k
B
. This is the Boltzmann distribution.
What is this mysterious
β
? Recall that the Lagrange multiplier
λ
E
measures
how S reacts to a change in
¯
E. In other words,
S
E
= λ
E
= k
B
β.
Moreover, by definition of temperature, we have
S
E
V,N,...
=
1
T
.
So it follows that
β =
1
k
B
T
.
Recall that the first law of thermodynamics says
dE = T dS P dV + µ dN + ··· .
This is a natural object to deal with when we have fixed
S, V, N
, etc. However,
often, it is temperature that is fixed, and it is more natural to consider the free
energy:
Definition
(Helmholtz free energy)
.
The Helmholtz free energy of a system at
fixed temperature, volume and particle number is defined by
F (T, V, N ) = U T S =
¯
E T S = k
B
T log Z.
This satisfies
dF = S dT P dV + µ dN + ··· ,
and is minimized at equilibrium for fixed T, V, N.
1.2 Coarse Graining
Usually, in statistical mechanics, we distinguish between two types of objects
microstates, namely the exact configuration of the system, and macrostates,
which are variables that describe the overall behaviour of the system, such that
pressure and temperature. Here we would like to consider something in between.
For example, if we have a system of magnets as in the Ising model, we the
microstate would be the magnetization at each site, and the macrostate would
be the overall magnetization. A coarse-graining of this would be a function
m
(
r
)
of space that describes the “average magnetization around
r
”. There is no fixed
prescription on how large an area we average over, and usually it does not matter
much.
In general, the coarse-grained variable would be called
ψ
. We can define a
coarse-grained partition function
Z[ψ(r)] =
X
iψ
e
βE
i
,
where we sum over all states that coarse-grain to
ψ
. We can similarly define the
energy and entropy by restricting to all such ψ, and get
F [ψ] = E[ψ] T S[ψ].
The probability of being in a state ψ is then
P[ψ] =
e
βF [ψ]
Z
TOT
, Z
TOT
=
Z
e
βF [ψ]
D[ψ].
What we have on the end is a functional integral, where we integrate over all
possible values of ψ. We shall go into details later. We then have
F
TOT
= k
B
T log Z
TOT
.
In theory, one can obtain
F
[
ψ
] by explicitly doing a coarse graining of the
macroscopic laws.
Example.
Consider an interacting gas with
N
particles. We can think of
the energy as a sum of two components, the ideal gas part (
d
2
NkT
), and an
interaction part, given by
E
int
=
1
2
X
i6j
U(r
i
r
j
),
where
i, j
range over all particles with positions
r
i
, r
j
respectively, and
U
is
some potential function. When we do coarse-graining, we introduce a function
ρ
that describes the local density of particles. The interaction energy can then be
written as
E
int
=
1
2
ZZ
U(r r
0
)ρ(r)ρ(r
0
) dr dr
0
.
Similarly, up to a constant, we can write the entropy as
S[ρ] = k
B
Z
ρ(r) log ρ(r) dr.
In practice, since the microscopic laws aren’t always accessible anyway, what
is more common is to take a phenomenological approach, namely we write down
a Taylor expansion of
F
[
ψ
], and then empirically figure out what the coefficients
should be, as a function of temperature and other parameters. In many cases, the
signs of the first few coefficients dictate the overall behaviour of the system, and
phase transition occurs when the change in temperature causes the coefficients
to switch signs.
2 Mean field theory
In this chapter, we explore the mean field theory of two physical systems
binary fluids and nematic liquid crystals. In mean field theory, what we do is we
write down the free energy of the system, and then find a state
φ
that minimizes
the free energy. By the Boltzmann distribution, this would be the “most likely
state” of the system, and we can pretend F
TOT
= F [φ].
This is actually not a very robust system, since it ignores all the fluctuations
about the minimum, but gives a good starting point for understanding the
system.
2.1 Binary fluids
Consider a binary fluid, consisting of a mixture of two fluids
A
and
B
. For
simplicity, we assume we are in the symmetric case, where
A
and
B
are the same
“type” of fluids. In other words, the potentials between the fluids are such that
U
AA
(r) = U
BB
(r) 6= U
AB
(r).
We consider the case where
A
and
B
repulse each other (or rather, repulse each
other more than the
A
-
A
and
B
-
B
repulsions). Thus, we expect that at high
temperatures, entropy dominates, and the two fluids are mixed together well. At
low temperatures, energy dominates, and the two fluids would be well-separated.
We let
ρ
A
(
r
) and
ρ
B
(
r
) be the coarse-grained particle density of each fluid,
and we set our order parameter to be
φ(r) =
ρ
A
(r) ρ
B
(r)
(N
A
+ N
B
)/V
,
with
N
A
and
N
B
, the total amount of fluids
A
and
B
, and
V
the volume. This
is normalized so that φ(r) [1, 1].
We model our system with Landau–Ginzburg theory, with free energy given
by
βF =
Z
a
2
φ
2
+
b
4
φ
4
| {z }
f(φ)
+
κ
2
(φ)
2
dr,
where a, b, κ are functions of temperature.
Why did we pick such a model? Symmetry suggests the free energy should
be even, and if we Taylor expand any even free energy functional, the first few
terms will be of this form. For small
φ
and certain values of
a, b, κ
, we shall see
there is no need to look further into higher order terms.
Observe that even without symmetry, we can always assume we do not have
a linear term, since a
term will integrate out to give
cV
¯
φ
, and
¯
φ
, the average
composition of the fluid, is a fixed number. So this just leads to a constant shift.
The role of the gradient term
R
κ
2
(
φ
)
2
d
r
captures at order
(2)
the non-
locality of E
int
,
E
int
=
X
i,j∈{A,B}
Z
ρ
i
(r)ρ
j
(r
0
)U
ij
(|r r
0
|) dr dr
0
,
If we assume
φ
(
r
) is slowly varying on the scale of interactions, then we can
Taylor expand this E
int
and obtain a (φ)
2
term.
Now what are the coefficients
a, b, κ
? For the model to make sense, we want
the free energy to be suppressed for large fluctuating
φ
. Thus, we want
b, κ >
0,
while
a
can take either sign. In general, the sign of
a
is what determines the
behaviour of the system, so for simplicity, we suppose
b
and
κ
are fixed, and let
a vary with temperature.
To do mean field theory, we find a single
φ
that minimizes
F
. Since the
gradient term
R
κ
2
(
φ
)
2
d
x
0, a naive guess would be that we should pick a
uniform φ,
φ(r) =
¯
φ.
Note that
¯
φ
is fixed by the constraint of the system, namely how much fluid of
each type we have. So we do not have any choice. In this configuration, the free
energy per unit volume is
F
V
= f(
¯
φ) =
a
2
¯
φ
2
+
b
4
¯
φ
4
.
The global of this function depends only on the sign of
a
. For
a >
0 and
a <
0
respectively, the plots look like this:
¯
φ
f
a > 0
a < 0
a > 0
a < 0
We first think about the a > 0 part. The key point is that the function f(φ) is
a convex function. Thus, for a fixed average value of
φ
, the way to minimize
f(φ) is to take φ to be constant. Thus, since
βF =
Z
f(φ(r)) +
κ
2
(φ)
2
dr,
even considering the first term alone tells us we must take
φ
to be constant, and
the gradient term reinforces this further.
The
a <
0 case is more interesting. The function
f
(
φ
) has two minima,
φ
1,2
= ±φ
B
, where
φ
B
=
r
a
b
.
Now suppose
¯
φ
lies between
±φ
B
. Then it might be advantageous to have some
parts of the fluid being at
φ
B
and the others at
φ
B
, and join them smoothly
in between to control the gradient term. Mathematically, this is due to the
concavity of the function f in the region [φ
B
, φ
B
].
Suppose there is
V
1
many fluid with
φ
=
φ
1
, and
V
2
many fluid with
φ
=
φ
2
.
Then these quantities must obey
V
1
φ
1
+ V
2
φ
2
= V
¯
φ,
V
1
+ V
2
= V.
Concavity tells us we must have
V
1
f(φ
1
) + V
2
f(φ
2
) < (V
1
+ V
2
)f(
¯
φ).
Thus, if we only consider the
f
part of the free energy, it is advantageous to
have this phase separated state. If our system is very large in size, since the
interface between the two regions is concentrated in a surface of finite thickness,
the gradient cost will be small compared to the gain due to phase separation.
We can be a bit more precise about the effects of the interface. In the first
example sheet, we will explicitly solve for the actual minimizer of the free energy
subject to the boundary condition
φ
(
x
)
±φ
B
as
x ±∞
, as in our above
scenario. We then find that the thickness of the interface is (of the order)
ξ
0
=
2κ
a
,
and the cost per unit area of this interface is
σ =
8κa
3
9b
2
1/2
.
This is known as the interfacial tension. When calculating the free energy of a
phase separated state, we will just multiply the interfacial tension by the area,
instead of going back to explicit free energy calculations.
In general the mean-field phase diagram looks like
a
1 1
¯
φ
a(T ) = 0
Within the solid lines, we have phase separation, where the ground state of the
system for the given
a
and
¯
φ
is given by the state described above. The inner
curve denotes spinodal instability, where we in fact have local instability, as
opposed to global instability. This is given by the condition
f
00
(
¯
φ
)
<
0, which
we solve to be
φ
S
=
r
a
3b
.
What happens if our fluid is no longer symmetric? In this case, we should
add odd terms as well. As we previously discussed, a linear term has no effect.
How about a cubic term
R
c
3
φ
(
r
)
3
d
r
to our
βF
? It turns out we can remove
the
φ
(
r
) term by a linear shift of
φ
and
a
, which is a simple algebraic maneuver.
So we have a shift of axes on the phase diagram, and nothing interesting really
happens.
2.2 Nematic liquid crystals
For our purposes, we can imagine liquid crystals as being made of rod-like
molecules
We are interested in the transition between two phases:
The isotropic phase, where the rods are pointed in random directions.
The nematic phase, where the rods all point in the same direction, so that
there is a long-range orientation order, but there is no long range positional
order.
In general, there can be two different sorts of liquid crystals the rods can
either be symmetric in both ends or have “direction”. Thus, in the first case,
rotating the rod by 180
does not change the configuration, and in the second
case, it does. We shall focus on the first case in this section.
The first problem we have to solve is to pick an order parameter. We want
to take the direction of the rod
n
, but mod it out by the relation
n n
. One
way to do so is to consider the second-rank traceless tensor
n
i
n
j
. This has the
property that
A
i
n
i
n
j
is the component of a vector
A
in the direction of
n
, and
is invariant under
n
i
7→ n
j
. Observe that if we normalize
n
to be a unit vector,
then
n
i
n
j
has trace 1. Thus, if we have isotropic rods in
d
dimensions, then we
have
hn
i
n
j
i =
δ
ij
d
.
In general, we can defined a coarse-grained order parameter to be
Q
ij
(r) = hn
i
n
j
i
local
1
d
δ
ij
.
This is then a traceless symmetric second-rank tensor that vanishes in the
isotropic phase.
One main difference from the case of the binary fluid is that
Q
ij
is no longer
conserved., i.e. the “total Q
Z
Q
ij
(r) dr
is not constant in time. This will have consequences for equilibrium statistical
mechanics, but also the dynamics.
We now want to construct the leading-order terms of the “most general” free
energy functional. We start with the local part
f
(
Q
), which has to be a scalar
built on Q. The possible terms are as follows:
(i) There is only one linear one, namely Q
ii
= Tr(Q), but this vanishes.
(ii)
We can construct a quadratic term
Q
ij
Q
ji
=
Tr
(
Q
2
), and this is in general
non-zero.
(iii)
There is a cubic term
Q
ij
Q
jk
Q
ki
=
Tr
(
Q
3
), and is also in general non-zero.
(iv) There are two possible quartic terms, namely Tr(Q
2
)
2
and Tr(Q
4
).
So we can write
f(Q) = a Tr(Q
2
) + c Tr(Q
3
) + b
1
Tr(Q
2
)
2
+ b
2
Tr(Q
4
).
This is the local part of the free energy up to fourth order in
Q
. We can go on,
and in certain conditions we have to, but if these coefficients
b
i
are sufficiently
positive in an appropriate sense, this is enough.
How can we think about this functional? Observe that if all of the rods point
tend to point in a fixed direction, say
z
, and are agnostic about the other two
directions, then Q will be given by
Q
ij
=
λ/2 0 0
0 λ/2 0
0 0 λ
, λ > 0.
If the rod is agnostic about the
x
and
y
directions, but instead avoids the
z
-direction, then
Q
ij
takes the same form but with
λ <
0. For the purposes of
f
(
Q
), we can locally diagonalize
Q
, and it should somewhat look like this form.
So this seemingly-special case is actually quite general.
The
λ >
0 and
λ <
0 cases are physically very different scenarios, but the
difference is only detected in the odd terms. Hence the cubic term is extremely
important here. To see this more explicitly, we compute f in terms of λ as
f(Q) = a
3
2
λ
2
+ c
3
4
λ
3
+ b
1
9
4
λ
4
+ b
2
9
8
λ
4
= ¯
2
+ ¯
3
+
¯
4
.
We can think of this in a way similar to the binary fluid, where
λ
is are sole
order parameter. We fix
¯
b
and
¯c <
0, and then vary
¯a
. In different situations,
we get
λ
f
α < α
c
α = α
c
α > α
c
Here the cubic term gives a discontinuous transition, which is a first-order
transition. If we had ¯c > 0 instead, then the minima are on the other side.
We now move on to the gradient terms. The possible gradient terms up to
order
(2)
and Q
(2)
are
κ
1
i
i
Q
j`
Q
j`
= κ
1
2
Tr(Q
2
)
κ
2
(
i
Q
im
)(
j
Q
jm
) = κ
2
( · Q)
2
κ
3
(
i
Q
jm
)(
j
Q
im
) = yuck.
Collectively, these three terms describe the energy costs of the following three
things:
splay
twisting
bend
In general, each of these modes correspond to linear combination of the three
terms, and it is difficult to pin down how exactly these correspondences work.
Assuming these linear combinations are sufficiently generic, a sensible choice is
to set
κ
1
=
κ
3
= 0 (for example), and then the elastic costs of these deformations
will all be comparable.
3 Functional derivatives and integrals
We shall be concerned with two objects functional derivatives and integrals.
Functional derivatives are (hopefully) familiar from variational calculus, and
functional integrals might be something new. They will be central to what we
are going to do next.
3.1 Functional derivatives
Consider a scalar field φ(r), and consider a functional
A[φ] =
Z
L(φ, φ) dr.
Under a small change
φ 7→ φ
+
δφ
(
r
) with
δφ
= 0 on the boundary, our functional
becomes
A[φ + δφ] =
Z
L(φ, φ) + δφ
L
φ
+ ·
L
φ
dr
= A[φ] +
Z
δφ
L
φ
·
L
φ
dr,
where we integrated by parts using the boundary condition. This suggests the
definition
δA
δφ(r)
=
L
φ(r)
·
L
φ
.
Example.
In classical mechanics, we replace
r
by the single variable
t
, and
φ
by position x. We then have
A =
Z
L(x, ˙x) dt.
Then we have
δA
δx(t)
=
L
x
d
dt
L
˙x
.
The equations of classical mechanics are
δA
δx(t)
= 0.
The example more relevant to us is perhaps Landau–Ginzburg theory:
Example. Consider a coarse-grained free energy
F [φ] =
Z
a
2
φ
2
+
b
4
φ
4
+
κ
2
(φ)
2
dr.
Then
δF
δφ(r)
= +
3
κ
2
φ.
In mean field theory, we set this to zero, since by definition, we are choosing
a single
φ
(
r
) that minimizes
F
. In the first example sheet, we find that the
minimum is given by
φ(x) = φ
B
tanh
x x
0
ξ
0
,
where ξ
0
is the interface thickness we previously described.
In general, we can think of
δF
δφ(r)
as a “generalized force”, telling us how we
should change
φ
to reduce the free energy, since for a small change
δφ
(
r
)), the
corresponding change in F is
δF =
Z
δF
δφ(r)
δφ(r) dr.
Compare this with the equation
dF = S dT p dV + µ dN + h · dM + ··· .
Under the analogy, we can think of
δF
δφ(r)
as the intensive variable, and
δφ
(
r
) as
the extensive variable. If
φ
is a conserved scalar density such as particle density,
then we usually write this as
µ(r) =
δF
δφ(r)
,
and call it the chemical potential. If instead
φ
is not conserved, e.g. the
Q
we
had before, then we write
H
ij
=
δF
δQ
ij
and call it the molecular field.
We will later see that in the case where
φ
is conserved,
φ
evolves according
to the equation
˙
φ = −∇ · J, J Dµ,
where
D
is the diffusivity. The non-conserved case is simpler, with equation of
motion given by.
˙
Q = ΓH.
Let us go back to the scalar field φ(r). Consider a small displacement
r 7→ r + u(r).
We take this to be incompressible, so that · u = 0. Then
φ 7→ φ
0
= φ
0
(r) = φ(r u).
Then
δφ(r) = φ
0
(r) φ(r) = u · φ(r) + O(u
2
).
Then
δF =
Z
δφ
δF
δφ
dr =
Z
µu · φ dr
=
Z
φ · (µu) dr =
Z
(φµ) · u dr =
Z
(φ
j
µ)u
j
dr.
using incompressibility.
We can think of the free energy change as the work done by stress,
δF =
Z
σ
ij
(r)ε
ij
(r) dr,
where
ε
ij
=
i
u
j
is the strain tensor, and
σ
ij
is the stress tensor. So we can
write this as
δF =
Z
σ
ij
i
u
j
dr =
Z
(
i
σ
ij
)u
j
dr.
So we can identify
i
σ
ij
= φ
j
µ.
So µ also contains the “mechanical information”.
3.2 Functional integrals
Given a coarse-grained ψ, we have can define the total partition function
e
βF
TOT
= Z
TOT
=
Z
e
βF [ψ]
D[ψ],
where D[
ψ
] is the “sum over all field configurations”. In mean field theory, we
approximate this
F
TOT
by replacing the functional integral by the value of the
integrand at its maximum, i.e. taking the minimum value of
F
[
ψ
]. What we
are going to do now is to evaluate the functional integral “honestly”, and this
amounts to taking into account fluctuations around the minimum (since those
far away from the minimum should contribute very little).
To make sense of the integral, we use the fact that the space of all
ψ
has a
countable orthonormal basis. We assume we work in [0
, L
]
q
of volume
V
=
L
q
with periodic boundary conditions. We can define the Fourier modes
ψ
q
=
1
V
Z
ψ(r)e
iq·r
dr,
Since we have periodic boundary conditions,
q
can only take on a set of discrete
values. Moreover, molecular physics or the nature of coarse-graining usually
implies there is some “maximum momentum”
q
max
, above which the wavelengths
are too short to make physical sense (e.g. vibrations in a lattice of atoms cannot
have wavelengths shorter than the lattice spacing). Thus, we assume
ψ
q
= 0 for
|q| > q
max
. This leaves us with finitely many degrees of freedom.
The normalization of ψ
q
is chosen so that Parseval’s theorem holds:
Z
|ψ|
2
dr =
X
q
|ψ
q
|
2
.
We can then define
D[ψ] =
Y
q
dψ
q
.
Since we imposed a
q
max
, this is a finite product of measures, and is well-defined.
In some sense,
q
max
is arbitrary, but for most cases, it doesn’t really matter
what
q
max
we choose. Roughly speaking, at really short wavelengths, the
behaviour of
ψ
no longer depends on what actually is going on in the system,
so these modes only give a constant shift to
F
, independent of interesting,
macroscopic properties of the system. Thus, we will mostly leave the cutoff
implicit, but it’s existence is important to keep our sums convergent.
It is often the case that after doing calculations, we end up with some
expression that sums over the
q
’s. In such cases, it is convenient to take the
limit V so that the sum becomes an integral, which is easier to evaluate.
An infinite product is still bad, but usually molecular physics or the nature
of coarse graining imposes a maximum
q
max
, and we take the product up to
there. In most of our calculations, we need such a
q
max
to make sense of our
integrals, and that will be left implicit. Most of the time, the results will be
independent of
q
max
(for example, it may give rise to a constant shift to
F
that
is independent of all the variables of interest).
Before we start computing, note that a significant notational annoyance is
that if
ψ
is a real variable, then
ψ
q
will still be complex in general, but they will
not be independent. Instead, we always have
ψ
q
= ψ
q
.
Thus, we should only multiply over half of the possible
q
’s, and we usually
denote this by something like
Q
+
q
.
In practice, there is only one path integral we are able to compute, namely
when βF is a quadratic form, i.e.
βF =
1
2
Z
φ(r)G(r r
0
)φ(r
0
) dr dr
0
Z
h(r)φ(r) dr.
Note that this expression is non-local, but has no gradient terms. We can think
of the gradient terms we’ve had as localizations of first-order approximations to
the non-local interactions. Taking the Fourier transform, we get
βF [ψ
q
] =
1
2
X
q
G(q)φ
q
φ
q
X
q
h
q
φ
q
.
Example. We take Landau–Ginzburg theory and consider terms of the form
βF [φ] =
Z
n
ξφ
2
+
κ
2
(φ)
2
+
γ
2
(
2
φ)
2
o
dr
The
γ
term is new, and is necessary because we will be interested in the case
where κ is negative.
We can now take the Fourier transform to get
βF {φ
q
} =
1
2
X
q
+
(a + κq
2
+ γq
4
)φ
q
φ
q
X
q
+
h
q
φ
q
.
=
X
q
+
(a + κq
2
+ γq
4
)φ
q
φ
q
X
q
+
h
q
φ
q
.
So our G(q) is given by
G(q) = a + kq
2
+ γq
4
.
To actually perform the functional integral, first note that if
h 6
= 0, then we
can complete the square so that the
h
term goes away. So we may assume
h
= 0.
We then have
Z
TOT
=
Z
"
Y
q
+
dφ
q
#
e
βF {φ
q
}
=
Y
q
+
Z
dφ
q
e
−|φ
q
|
2
G(q)
Each individual integral can be evaluated as
Z
dφ
q
e
−|φ
q
|
2
G(q)
=
Z
ρ dρ dθ e
G(q)ρ
2
=
π
G(q)
,
where φ
q
= ρe
. So we find that
Z
TOT
=
Y
q
+
π
G(q)
,
and so
βF
T
= log Z
T
=
X
q
+
log
G(q)
π
.
We now take large V limit, and replace the sum of the integral. Then we get
βF
T
=
1
2
V
(2π)
d
Z
q
max
dq log
G(q)
π
.
There are many quantities we can compute from the free energy.
Example. The structure factor is defined to be
S(k) = hφ
k
φ
k
i =
1
Z
T
Z
φ
k
φ
k
e
P
+
q
φ
q
φ
q
G(q)
Y
q
+
dφ
q
.
We see that this is equal to
1
Z
T
Z
T
G(k)
=
log Z
T
G(k)
=
1
G(k)
.
We could also have done this explicitly using the product expansion.
This
S
(
k
) is measured in scattering experiments. In our previous example,
for small k and κ > 0, we have
S(q) =
1
a + κk
2
+ γk
4
a
1
1 + k
2
ξ
2
, ξ =
r
κ
a
.
where ξ is the correlation length. We can return to real space by
hφ
2
(r)i =
1
V
Z
|φ(r)|
2
dr
=
1
V
X
q
hφ
q
|
2
i
=
1
(2π)
d
Z
q
max
dq
a + κq
2
+ γq
4
.
4 The variational method
4.1 The variational method
The variational method is a method to estimate the partition function
e
βF
TOT
=
Z
e
βF [φ]
D[φ]
when
F
is not Gaussian. To simplify notation, we will set
β
= 1. It is common
to make a notational change, where we replace
F
TOT
with
F
and
F
[
φ
] with
H
[
φ
].
We then want to estimate
e
F
TOT
=
Z
e
F [φ]
D[φ].
We now make a notation change, where we write
F
TOT
as
F
, and
F
[
φ
] as
H
[
φ
]
instead, called the effective Hamiltonian. In this notation, we write
e
F
=
Z
e
H[φ]
D[φ].
The idea of the variational method is to find some upper bounds on
F
in
terms of path integrals we can do, and then take the best upper bound as our
approximation to F .
Thus, we introduce a trial Hamiltonian H
0
[φ], and similarly define
e
F
0
=
Z
e
H
0
[φ]
D[φ].
We can then write
e
F
=
e
F
0
R
e
H
0
D[φ]
Z
e
H
0
e
(HH
0
)
D[φ] = e
F
0
he
(HH
0
)
i
0
,
where the subscript 0 denotes the average over the trial distribution. Taking the
logarithm, we end up with
F = F
0
loghe
(HH
0
)
i
0
.
So far, everything is exact. It would be nice if we can move the logarithm inside
the expectation to cancel out the exponential. While the result won’t be exactly
equal, the fact that log is concave, i.e.
log(αA + (1 α)B) α log A + (1 α) log B.
Thus Jensen’s inequality tells us
loghY i
0
hlog Y i
0
.
Applying this to our situation gives us an inequality
F F
0
hH
0
Hi
0
= F
0
hH
0
i
0
+ hHi
0
= S
0
+ hHi
0
.
This is the Feynman–Bogoliubov inequality.
To use this, we have to choose the trial distribution
H
0
simple enough to
actually do calculations (i.e. Gaussian), but include variational parameters in
H
0
. We then minimize the quantity
F
0
hH
0
i
+
hHi
0
over our variational
parameters, and this gives us an upper bound on
F
. We then take this to be
our best estimate of
F
. If we are brave, we can take this minimizing
H
0
as an
approximation of H, at least for some purposes.
4.2 Smectic liquid crystals
We use this to talk about the isotropic to smectic transition in liquid crystals.
The molecules involved often have two distinct segments. For example, we may
have soap molecules that look like this:
The key property of soap molecules is that the tail hates water while the head
likes water. So we expect these molecules to group together like
In general, we can imagine our molecules look like
and like attracts like. As in the binary fluid, we shall assume the two heads are
symmetric, so
U
AA
=
U
BB
6
=
U
AB
. If we simply want the different parts to stay
away from each other, then we can have a configuration that looks like
z
In general, we expect that there is such an order along the
z
direction, as
indicated, while there is no restriction on the alignments in the other directions.
So the system is a lattice in the
z
direction, and a fluid in the remaining two
directions. This is known as a smectic liquid crystal, and is also known as the
lamellar phase. This is an example of microphase separation.
As before, we let
φ
be a coarse grained relative density. The above ordered
phase would then look like
φ(x) = cos q
0
z.
for some
q
0
that comes from the molecular length. If our system is not perfectly
ordered, then we may expect it to look roughly like A cos q
0
z for some A.
We again use Landau–Ginzburg model, which, in our old notation, has
βF =
Z
a
2
φ
2
+
b
4
φ
4
+
κ
2
(φ)
2
+
γ
2
(
2
φ)
2
dr.
If we write this in Fourier space, then we get
βF =
1
2
X
q
(a + κq
2
+ γq
4
)φ
q
φ
q
+
b
4V
X
q
1
,q
2
,q
3
φ
q
1
φ
q
2
φ
q
3
φ
(q
1
+q
2
+q
3
)
.
Notice that the quartic term results in the rather messy sum at the end. For the
iso-smectic transition, we choose κ < 0, γ > 0.
Again for simplicity, we first consider the case where
b
= 0. Then this is a
Gaussian model with
G(q) = a + κq
2
+ γq
4
= τ + α(q q
0
)
2
,
Varying
a
gives a linear shift in
G
(
q
). As we change
a
, we get multiple different
curves.
q
G(q)
a = a
c
q
0
Thus, as a decreases, S(q) = h|φ
q
|
2
i =
1
G(q)
blows up at some finite
q = q
0
=
r
κ
2γ
, a
c
=
κ
2
4γ
.
We should take this blowing up as saying that the
|q|
=
q
0
states are highly
desired, and this results in an ordered phase. Note that any
q
with
|q|
=
q
0
is
highly desired. When the system actually settles to an ordered state, it has to
pick one such
q
and let the system align in that direction. This is spontaneous
symmetry breaking.
It is convenient to complete to square, and expand
G
about
q
=
q
0
and
a = a
c
. Then we have
G(q) = τ + α(q q
0
)
2
,
where
τ = a a
c
, α =
1
2
G
00
(q
0
) = 2κ.
Then the transition we saw above happens when τ = 0.
We now put back the quartic term. We first do this with mean field theory,
and then later return to the variational method.
Mean Field Theory
In mean field theory, it is easier to work in real space. We look for a single field
configuration that minimizes
F
. As suggested before, we try a solution of the
form
φ = A cos q
0
z,
which is smectic along
z
. We can then evaluate the free energy (per unit volume)
to be
βF
V
= βF [φ] =
a
2
φ
2
+
κ
2
(φ)
2
+
γ
2
(
2
φ)
2
+
b
4
φ
4
where the bar means we average over one period of the periodic structure. It is
an exercise to directly compute
φ
2
=
1
2
A
2
, (φ)
2
=
1
2
A
2
q
2
0
, (
2
φ)
2
=
1
2
A
2
q
4
0
, φ
4
=
3
8
A
4
.
This gives
βF
V
=
1
4
aA
2
+ κA
2
q
2
0
+ γA
2
q
4
0
+
3b
8
A
4
=
1
4
A
2
(a a
c
)
| {z }
τ
+
3b
8
A
4
.
Note that
q
0
is fixed by the system as we originally had, while
A
is the amplitude
of the fluctuation which we get to choose. Plotting this, we get a familiar graph
A
βF/V
τ > 0
τ < 0
τ > 0
τ < 0
If
τ >
0, then the optimum solution is given by
A
= 0. Otherwise, we should
pick
A 6
= 0. Observe that as we slowly reduce
τ
across 0, the minimum varies
continuously with τ :
a
|A|
a
c
Variational method
We now consider the variational theory. In the notation we were using previously,
we have
H =
X
q
+
φ
q
φ
q
G(q) +
b
4V
X
q
φ
q
1
φ
q
2
φ
q
3
φ
(q
1
+q
2
+q
3
)
.
Our trial H
0
is
H
0
=
X
q
φ
q
φ
q
J(q).
Since this is a Gaussian model, we know that
F
0
=
X
q
+
log
J(q)
π
.
To use our inequality, we need to evaluate our other two bits. We have
hH
0
i
0
=
X
q
+
hφ
q
φ
q
i
0
J(q).
We already calculated
hφ
q
φ
q
i
0
=
1
J(q)
.
Thus, we have
hH
0
i
0
=
X
q
+
1.
Here it is clear that we must impose a cutoff on
q
. We can think of this 1 as the
equipartition theorem.
We can also compute
hHi
0
=
X
q
+
1
J(q)
G(q) +
b
4V
X
q
1
,q
2
,q
3
hφ
q
1
φ
q
2
φ
q
3
φ
(q
1
q
2
q
3
)
i
0
| {z }
.
In the Gaussian model, each
φ
q
is a zero mean Gaussian random variables, which
have certain nice properties. Wick’s theorem tells us we have
habcdi
0
= habi
0
hcdi
0
+ haci
0
hbdi
0
+ hadi
0
hbci
0
.
Applying this, together with the result
hφ
q
1
φ
q
0
i
0
= h|φ
q
|
2
i
0
δ
q
1
,q
2
,
we obtain
U = 3
"
X
q
h|φ
q
|
2
i
0
#
2
= 12
"
X
q
+
1
J(q)
#
2
.
Thus, we have
hH
0
i
0
=
X
q
+
G(q)
J(q)
+
3b
V
X
q
+
1
J(q)
!
2
.
We can then compute
˜
F =
X
q
+
log
J(q)
π
1 +
G(q)
J(q)
+
3b
V
X
q
+
1
J(q)
!
2
.
We minimize over J(q) by solving
˜
F
J(q)
= 0
for all q. Differentiating, we obtain
1
J(q)
G(q)
J(q)
2
6b
V J(q)
2
X
q
0
+
1
J(q
0
)
= 0.
Multiplying through by J
2
, we get
J(q) = G(q) +
6b
V
X
q
0
+
1
J(q
0
)
.
For large V , we can replace the sum by an integral, and we have
J(q) = G(q) +
3b
(2π)
d
Z
dq
0
J(q
0
)
.
It is very important that once we have fixed
J
(
q
), the second term is a constant.
Writing
C =
2
V
X
q
0
+
1
J(q
0
)
=
1
(2π)
d
Z
dq
0
J(q
0
)
,
we can then say the minimum is given by
J
(
q
) =
G
(
q
) + 3
bC
. Thus, solving for
J(q) is equivalent to finding a C such that
C =
1
(2π)
d
Z
dq
0
G(q) + 3bC
.
This is a self-consistency equation for C (and hence J).
There is a very concrete interpretation of
C
. Recall that
hφ
2
i
is the average
value of φ
2
at some point r (which is independent of r. So we can compute
hφ(r)
2
i =
1
V
Z
hφ(r)
2
i dr =
1
V
X
q
h|φ
q
|
2
i =
1
V
X
q
1
J(q)
.
Thus, what our computation above amounts to is that we replaced
G
(
q
) by
G(q) + 3bhφ
2
i.
A good way to think about this is that we are approximating the
φ
4
term in
the free energy by
b
4
φ
4
3b
2
hφ
2
iφ
2
.
We have to pick a value so that this approximation is consistent. We view the
hφ
2
i
above as just another constant, so that the free energy is now a Gaussian.
We can compute the expectation of
φ
2
using this free energy, and then we find
that
hφ
2
i =
1
(2π)
d
Z
dq
0
G(q) + 3bhφ
2
i
.
This is the self-consistency equation as above.
We could have done the approximation above without having to through the
variational method, but the factor of
3b
2
is then no longer obvious.
To solve this self-consistency equation, or at least understand it, it is conve-
nient to think in terms of τ instead. If we write our G(q) as
G(q) = τ + α(q q
0
)
2
,
then we have
J(q) = ¯τ + α(q q
0
)
2
, ¯τ = 3bhφ
2
i.
The self-consistency equation now states
¯τ = τ +
3b
(2π)
d
Z
d
d
q
¯τ + α(q q
0
)
2
.
We are interested in what happens near the critical point. In this case, we expect
¯τ
to be small, and hence the integral is highly peaked near 0. In
d
= 3, we can
make the approximation.
3b
(2π)
3
Z
d
3
q
¯τ + α(q q
0
)
2
=
3b
2π
2
Z
0
q
2
dq
¯τ + α(q q
0
)
2
3bq
2
0
2π
2
Z
0
dq
¯τ + α(q q
0
)
2
.
While there is a nasty integral to evaluate, we can make the substitution
q 7→ ¯τq
to bring the dependency on ¯τ outside the integral, and obtain
¯τ = τ +
sb
¯τ
, s =
3q
2
0
2π
2
1
α
Z
0
dy
1 + y
2
q
2
0
α
.
The precise value of
s
does not matter. The point is that it is constant indepen-
dent of
¯τ
. From this, we see that
¯τ
can never reach zero! This means we can
never have a continuous transition. At small
¯τ
, we will have large fluctuations,
and the quantity
h|φ
q
|
2
i = S(q) =
1
¯τ + α(q q
0
)
2
becomes large but finite. Note that since
¯τ
is finite, this sets some sort of length
scale where all q with |q q
0
| ¯τ have large amplitudes.
We can think of the above computation as looking at neighbourhoods of the
φ
= 0 vacuum, as that is what the variational methods see. So
¯τ
never vanishes
in this regime, we would expect a discontinuous isotropic-smectic transition that
happens “far away” from this vacuum.
Consider fluctuations about an ordered state, φ = φ
0
+ δφ, where
φ
0
= A cos q
0
z.
We can do computations similar to what we did above by varying
δφ
, and obtain
a
˜
F
(
A
). Then the global minimum over all
A
then gives us the “true” ground
state. Instead of doing that, which is quite messy, we use the heuristic version
instead. For A = 0, we had
τ = ¯τ + 3bhφ(r)
2
i.
For finite A, the quartic term now has an extra contribution, and we bet
¯τ = τ +
sb
¯τ
+
3bA
2
2
.
Compare this with mean field theory, where we have
F
MF T
(A) =
V
4
τA
2
+
3b
8
A
4
.
We see that for small
A
, the fluctuations are large, and mean field theory is
quite off. For large
A
, the fluctuation terms are irrelevant and mean field theory
is a good approximation. Thus, we get
A
F
MFT
We see that as
τ
decreases, the minimum discontinuously jumps from
A
= 0
to a finite value of A. Mean field theory is approached at large A.
We can then plot the minimum value of A as a function of τ:
MFT
τ
|A|
τ
c
A
c
τ
H
We have not calculated
A
c
or
τ
c
, but it can be done. We shall not do this, but
we can state the result (Brazovskii, 1975):
τ
c
' (sb)
2/3
, A
c
' s
1/3
b
1/6
.
It turns out the variational approach finally breaks down for
τ < τ
H
s
3/4
b
1/2
.
We have
τ
H
τ
c
if
b
s
. The reason this breaks down is that at low enough
temperatures, the fluctuations from the quartic terms become significant, and
our Gaussian approximation falls apart.
To figure out what
τ
H
is, we need to find leading corrections to
˜
F
, as
Brazovskii did. In general, there is no reason to believe
τ
H
is large or small. For
example, this method breaks down completely for the Ising model, and is correct
in no regime at all. In general, the self-consistent approach here is ad hoc, and
one has to do some explicit error analysis to see if it is actually good.
Brazovskii transition with cubic term
What happens when we have a cubic term? This would give an
A
3
term at
mean field level, which gives a discontinuous transition, but in a completely
different way. We shall just state the results here, using a phase diagram with
two parameters τ and c. In mean field theory, we get
c
τ
I
S
H
where H is a hexagonal phase, which looks like
where each cylinder contains a high concentration of a fixed end of the molecule.
This is another liquid crystal, with two crystal directions and one fluid directions.
The self-consistent version instead looks like
c
τ
I
S
H
˜c
Here c only matters above a threshold ˜c.
5 Dynamics
We now want to understand dynamics, namely if we have a system out of
equilibrium, how will it evolve in time? Physically, such situations are usually
achieved by rapidly modifying external parameters of a system. For example,
if the system is temperature-dependent, one may prepare a sample at high
temperature so that the system is in a homogeneous state, and then quench the
system by submerging it in water to lower the temperature rapidly. The system
will then slowly evolve towards equilibrium.
Before we think about the problem of dynamics, let’s think about a more
fundamental question what is it that is preventing the system from collapsing
to the ground state entirely, as opposed to staying in the Boltzmann distribution?
The answer is that our system is in contact with a heat bath, which we can
model as some random noise driving the movement of our particles. This gives a
dynamical way of achieving the Boltzmann distribution. When the system is
out of equilibrium, the random noise is still present and drives our system. The
key point is that the properties of the noise can be derived from the fact that at
equilibrium, they give the Boltzmann distribution.
5.1 A single particle
We ultimately want to think about field theories, but it is helpful to first consider
the case of a single, 1-dimensional particle. The action of the particle is given by
A =
Z
L dt, L = T V.
The equations of motion are
δA
δx(t)
=
L
x
d
dt
L
x
= 0.
For example, if
L =
1
2
m ˙x
2
V (x),
then the equation of motion is
δA
δx(t)
= m¨x + V
0
(x) = 0.
There are two key properties of this system
The system is deterministic.
The Lagrangian is invariant under time reversal. This is a consequence of
the time reversal symmetry of microscopic laws.
We now do something different. We immerse the particle in a fluid bath,
modelling the situation of a colloid. If we were an honest physicist, we would
add new degrees of freedom for each of the individual fluid particles. However,
we are dishonest, and instead we aggregate these effects as new forcing terms in
the equation of motion:
(i) We introduce damping, F
D
= ζ ˙x.
(ii) We introduce a noise f, with hfi = 0.
We set F
BATH
= F
D
+ f . Then we set our equation of motion to be
δA
δx
= F
BATH
ζ ˙x + f.
This is the Langevin equation.
What more can we say about the noise? Physically, we expect it to be the
sum of many independent contributions from the fluid particles. So it makes
sense to assume it takes a Gaussian distribution. So the probability density of
the realization of the noise being f is
P[f(t)] = N
f
exp
1
σ
2
Z
f(t)
2
dt
,
where
N
f
is a normalization constant and
σ
is the variance, to be determined.
This is called white noise. This has the strong independence property that
hf(t)f(t
0
)i = σ
2
δ(t t
0
).
In this course, we always assume we have a Gaussian white noise.
Since we have a random noise, in theory, any path we can write down is a
possible actual trajectory, but some are more likely than others. To compute
the probability density, we fixed start and end points (
x
1
, t
1
) and (
x
2
, t
2
). Given
any path
x
(
t
) between these points, we can find the noise of the trajectory to be
f = ζ ˙x
δA
δx
.
We then can compute the probability of this trajectory happening as
P
F
[x(t)] “=” P[f] = N
x
exp
1
2σ
2
Z
t
2
t
1
ζ ˙x
δA
δx
2
dt
!
.
This is slightly dodgy, since there might be some Jacobian factors we have missed
out, but it doesn’t matter at the end.
We now consider the problem of finding the value of
σ
. In the probability
above, we wrote it as
P
F
, denoting the forward probability. We can also consider
the backward probability
P
B
[
x
(
t
)], which is the probability of going along the
path in the opposite direction, from (x
2
, t
2
) to (x
1
, t
1
).
To calculate this, we use the assumption that
δA
δx
is time-reversal invariant,
whereas ˙x changes sign. So the backwards probability is
P
B
[x(t)] = N
x
exp
1
2σ
2
Z
t
2
t
1
ζ ˙x
δA
δx
2
!
dt.
The point is that at equilibrium, the probability of seeing a particle travelling
along
x
(
t
) forwards should be the same as the probability of seeing a particle
travelling along the same path backwards, since that is what equilibrium means.
This is not the same as saying
P
B
[
x
(
t
)] is the same as
P
F
[
x
(
t
)]. Instead, if at
equilibrium, there are a lot of particles at
x
2
, then it should be much less likely
for a particle to go from x
2
to x
1
than the other way round.
Thus, what we require is that
P
F
[x(t)]e
βH
1
= P
B
[x(t)]e
βH
2
,
where
H
i
=
H
(
x
i
, ˙x
i
). This is the principle of detailed balance. This is a
fundamental consequence of microscopic reversibility. It is a symmetry, so coarse
graining must respect it.
To see what this condition entails, we calculate
P
F
[x(t)]
P
B
[x(t)]
=
exp
1
2σ
2
R
t
2
t
1
(ζ ˙x)
2
2ζ ˙x
δA
δx
+
δA
δx
2
dt
exp
1
2σ
2
R
t
2
t
1
(ζ ˙x)
2
+ 2ζ ˙x
δA
δx
+
δA
δx
2
dt
= exp
1
2σ
2
Z
t
2
t
1
4ζ ˙x
δA
δx
dt
.
To understand this integral, recall that the Hamiltonian is given by
H(x, ˙x) = ˙x
L
˙x
L.
In our example, we can explicitly calculate this as
H =
1
2
˙x
2
+ V (x).
We then find that
dH
dt
=
d
dt
˙x
L
˙x
dL
dt
= ¨x
L
˙x
+ ˙x
d
dt
L
˙x
˙x
L
x
+ ¨x
L
˙x
= ˙x
d
dt
L
˙x
L
x
= ˙x
δA
δx
.
Therefore we get
Z
t
2
t
1
˙x
δA
δx(t)
dt = (H(x
2
, ˙x
2
) H(x
1
, ˙x
1
)) = (H
2
H
1
).
Therefore, the principle of detailed balance tells us we should pick
σ
2
= 2k
B
T ζ.
This is the simplest instance of the fluctuation dissipation theorem.
Given this, we usually write
f =
p
2k
B
T ζΛ,
where Λ is a Gaussian process and
hΛ(t)Λ(t
0
)i = δ(t t
0
).
We call Λ a unit white noise.
In summary, for a particle under a potential V , we have an equation
m¨x + V
0
(x) = ζ ˙x + f.
The term
ζ ˙x
gives an arrow of time en route to equilibrium, while the noise
term resolves time reversal symmetry once equilibrium is reached. Requiring
this fixes the variance, and we have
hf(t)f(t
0
)i = σ
2
δ(t t
0
) = 2k
B
T ζδ(t t
0
).
In general, in the coarse grained world, suppose we have mesostates
A, B
,
with probabilities e
βF
A
and e
βF
B
, then we have an identical statement
e
βF
A
P(A B) = e
βF
B
P(B A).
5.2 The Fokker–Planck equation
So far, we have considered a single particle, and considered how it evolved over
time. We then asked how likely certain trajectories are. An alternative question
we can ask is if we have a probability density for the position of
x
, how does
this evolve over time?
It is convenient to consider the overdamped limit, where
m
= 0. Our equation
then becomes
ζ ˙x = −∇V +
p
2k
B
T ζΛ.
Dividing by ζ and setting
˜
M = ζ
1
, we get
˙x =
˜
MV +
q
2k
B
T
˜
MΛ.
This
˜
M is the mobility, which is the velocity per unit force.
We define the probability density function
P (x, t) = probability density at x at time t.
We can look at the probability of moving by a distance
x
in a time interval
t. Equivalently, we are asking Λ to change by
∆Λ =
1
2k
B
T ζ
(ζx + V · t).
Thus, the probability of this happening is
W (∆x, x) P
t
(∆x) = N exp
1
4ζk
B
T t
(ζx + V t)
2
.
We will write
u
for
x
. Note that
W
(
u, x
) is just a normal, finite-dimensional
Gaussian distribution in
u
. We can then calculate that after time
t
, the
expectation and variance of u are
hui =
V
ζ
t, hu
2
i hui
2
=
2k
B
T
ζ
t + O(∆t
2
).
We can find a deterministic equation for P (x, t), given in integral form by
P (x, t + t) =
Z
P (x u, t)W (u, x u) du.
To obtain a differential form, we Taylor expand the integrand as
P (x u, t)W (u, x u)
=
P uP +
1
2
u
2
2
P
W (u, x) uW +
1
2
u
2
2
W
,
where all the gradients act on
x
, not
u
. Applying the integration to the expanded
equation, we get
P (x, t + t) = P (x, t) hui∇P +
1
2
hu
2
i∇
2
P P ∇hui,
Substituting in our computations for hui and hu
2
i gives
˙
P (x, t)∆t =
V
ζ
P +
k
B
T
ζ
2
P +
1
ζ
P
2
V
t.
Dividing by t, we get
˙
P =
k
B
T
ζ
2
P +
1
ζ
(P V )
= D
2
P +
˜
M(P V ),
where
D =
k
B
T
ζ
,
˜
M =
D
k
B
T
=
1
ζ
are the diffusivity and the mobility respectively.
Putting a subscript
1
to emphasize that we are working with one particle,
the structure of this is
˙
P
1
= −∇ · J
1
J
1
= P
1
D(log P
1
+ βV )
= P
1
˜
Mµ(x),
where
µ = k
B
T log P
1
+ V
is the chemical potential of a particle in V (x), as promised. Observe that
This is deterministic for P
1
.
This has the same information as the Langevin equation, which gives the
statistics for paths x(t).
This was “derived” for a constant
ζ
, independent of position. However,
the equations in the final form turn out to be correct even for
ζ
=
ζ
(
x
) as
long as the temperature is constant, i.e.
˜
M
=
˜
M
(
x
) =
D(x)
k
B
T
. In this case,
the Langevin equation says
˙x =
˜
M(x)V +
p
2D(x.
The multiplicative (i.e. non-constant) noise term is problematic. To
understand multiplicative noise, we need advanced stochastic calculus
(Itˆo/Stratonovich). In this course (and in many research papers), we avoid
multiplicative noise.
5.3 Field theories
Suppose now that we have
N
non-interacting colloids under the same potential
V
(
x
). As usual, we model this by some coarse-grained density field
ρ
(
x, t
). If
we assume that the particles do not interact, then the expected value of this
density is just given by
hρi = N P
1
, h˙ρi = N
˙
P
1
.
Then our previous discussion entails ˙ρ evolves by
˙ρ = −∇ · J,
where
hJi = ρ
˜
Mµ, µ = k
B
T log ρ + V (x).
If we wish to consider a general, interacting field, then we can take the same
equations, but set
µ =
δF
δρ
instead.
Note that these are hydrodynamic level equations for
ρ
, i.e. they only tell
us what happens to
hρi
. If we put
J
=
hJi
, then we get a mean field solution
that evolves to the minimum of the free energy. To understand the stochastic
evolution of ρ itself, we put
J = ρ
˜
Mµ + j,
where
j
is a noise current. This is the Langevin equation for a fluctuating field
ρ(r, t).
We can fix the distribution of
j
by requiring detailed balance as before. We
will implement this for a constant
M
=
ρ
˜
M
, called the collective mobility. This
is what we have to do to avoid having multiplicative noise in our system. While
this doesn’t seem very physical, this is reasonable in situations where we are
looking at small fluctuations about a fixed density, for example.
As before, we assume j(r, t) is Gaussian white noise, so
P[j(r, t)] = N exp
1
2σ
2
Z
t
2
t
1
dt
Z
dr |j(r, t)|
2
.
This corresponds to
hj
k
(r, t)j
`
(r
0
, t
0
)i = σ
2
δ
k`
δ(r r
0
)δ(t t
0
).
We now repeat the detailed balance argument to find σ
2
. We start from
J + Mµ = j.
Using F to mean forward path, we have
P
F
[J(r, t)] = N exp
1
2σ
2
Z
t
2
t
1
dt
Z
dr |J + Mµ|
2
,
where
µ =
δF [ρ]
δρ
.
We consider the backwards part and get
P
B
[J(r, t)] = N exp
1
2σ
2
Z
t
2
t
1
dt
Z
dr | J + M µ|
2
,
Then
log
P
F
P
B
=
2M
σ
2
Z
t
2
t
1
dt
Z
dr J · µ.
We integrate by parts in space to write
Z
dr J · µ =
Z
dr ( · J)µ =
Z
dr
˙ρ
δF
δρ
=
dF [ρ]
dt
.
So we get
log
P
F
P
B
=
2M
σ
2
(F
2
F
1
).
So we need
2M
σ
2
= β,
or equivalently
σ
2
= 2k
B
T M.
So our final many-body Langevin equation is
˙ρ = −∇ · J
J = M
δF
δρ
+
p
2k
B
T MΛ,
where Λ is spatiotemporal unit white noise. As previously mentioned, a constant
M avoids multiplicative white noise.
In general, we get the same structure for any other diffusive system, such as
φ(r, t) in a binary fluid.
We might want to get a Fokker–Planck equation for our field theory. First
recap what we did. For one particle, we had the Langevin equation
˙x =
˜
MV +
q
2k
B
T
˜
MΛ,
and we turned this into a Fokker–Planck equation
˙
P = −∇· J
J = P
˜
Mµ
µ = k
B
T log P + V (x).
We then write this as
˙
P = ·
˜
Mk
B
T ( + βV ) P
,
where P (x, t) is the time dependent probability density for x.
A similar equation can be derived for the multi-particle case, which we will
write down but not derive. We replace
x
(
t
) with
ρ
(
r, t
), and we replace
P
(
x, t
)
with
P
[
ρ
(
r
);
t
]. We then replace
with
δ
δρ(r)
. So the Fokker–Planck equation
becomes
˙
P [ρ(t); t] =
Z
dr
δ
δρ
k
B
T ·
˜
M
δ
δρ
+ β
δF
δρ
P
.
This is the Fokker–Planck equation for fields ρ.
As one can imagine, it is not very easy to solve. Note that in both cases, the
quantities
+
βV
and
δ
δρ
+
β
δF
δρ
annihilate the Boltzmann distribution. So
the Boltzmann distribution is invariant.
The advantage of the Langevin equation is that it is easy to understand the
mean field theory/deterministic limit
ρ
=
ρ
hydro
(
r, t
). However, it is difficult to
work with multiplicative noise. In the Fokker–Planck equation, multiplicative
noise is okay, but the deterministic limit may be singular. Schematically, we
have
P [ρ(r), t] = δ(ρ(r, t) ρ
hydro
(r, t)).
In this course, we take the compromise and use the Langevin equation with
constant M.
6 Model B
We now apply this theory to model some concrete systems. We shall first consider
a simple model of binary fluids. We assume that diffusion happens but without
fluid flow. As before, this is modeled by a scalar composition field
φ
(
r
), and
evolves under the equations
˙
φ = −∇ · J
J = M µ +
p
2k
B
T MΛ
µ =
δF
δφ
.
Recall that the system is modelled under the Landau–Ginzburg free energy
F [φ] =
Z
a
2
φ
2
+
b
4
φ
4
| {z }
f
+
κ
2
(φ)
2
dr.
We then have
µ = +
3
κ
2
φ.
As before, the mean field theory for equilibrium looks like
¯
φ
f
a > 0
a < 0
a > 0
a < 0
φ
B
φ
S
Here
±φ
S
are the spinodals, where the second derivative changes sign. This
gives the phase diagram
a
1 1
¯
φ
a(T ) = 0
1
3
2
Here
¯
φ is the global composition, which is a control variable.
In the past, we discussed what the field looks like in each region when we are
at equilibrium. At (1), the system is in a uniform phase that is globally stable.
If we set up our system at (1), and then rapidly change the temperature so that
we lie in (2) or (3), then we know that after the system settles, we will have a
phase separated state. However, how this transition happens is not something
mean field theory tells us. Heuristically, we expect
In (2), we have
|
¯
φ| < φ
S
, and
f
00
(
φ
)
<
0 implies local instability. The
system rapidly becomes phase separated. This is spinodal behaviour.
In (3), we have
φ
S
< |
¯
φ| < φ
B
. A uniform phase is locally stable, but
not globally. To reach the phase separated state, we need nucleation and
growth to occur, and requires the contribution of noise.
We now study these in detail.
Regime 1
We know that regime (1) is stable, and we shall see how it responds to perturba-
tion about φ(r) =
¯
φ. Put
φ =
¯
φ +
˜
φ(r).
We can then write
µ =
δF
δφ
=
f
φ
κ
2
φ = f
0
(
¯
φ) +
˜
φf
00
(
¯
φ) κ
2
˜
φ.
Note that the first term is a constant. We then have
˙
φ = −∇ · J
J = M [f
00
˜
φ κ
2
˜
φ] +
p
2k
B
T MΛ.
We drop the tildes and take the Fourier transform to get
˙
φ
q
= Mq
2
(f
00
+ κq
2
)φ
q
+ iq ·
p
2k
B
T MΛ
q
.
Compare this with an overdamped particle in a simple harmonic oscillator,
V =
1
2
κx
2
,
where we have
˙x =
˜
Mκx +
q
2k
B
T
˜
MΛ.
Indeed, we can think of our system as an infinite family of decoupled harmonic
oscillators, and solve each of them independently.
In the second example sheet, we compute
S(q, t) hφ
q
(0)φ
q
(t)i = S(q)e
r(q)t
,
This
S
(
q, t
) is called the dynamic structure factor , which can be measured by light
scattering. This doesn’t say the fluctuations go away completely we expect
there to be fluctuations all the time. What this says is that fluctuations at late
times come completely from the random noise, and not the initial fluctuations.
Regime 2
Consider the case where we are in the second regime. As before, we have the
equation
˙
φ
q
= Mq
2
(f
00
+ κq
2
)
| {z }
r(q)
φ
q
+ iq ·
p
2k
B
T MΛ
q
,
but crucially, now
f
00
(
¯
φ
)
<
0, so it is possible to have
r
(
q
)
<
0. The system is
unstable.
If we ignore the noise by averaging the equation, then we get
h
˙
φ
q
i = r(q)hφ
q
i.
So if we have a noisy initial state φ
q
(0), then the perturbation grows as
hφ
q
(t)i = φ
q
(0)e
r(q)t
.
When
r
(
q
)
<
0, then this amplifies the initial noise. In this world, even if we
start with a perfectly uniform
φ
, noise terms will kick in and get amplified over
time. Moreover, since we have an exponential growth, the earliest noise gets
amplified the most, and at late times, the new perturbations due to the noise
are negligible.
We can plot our r(q) as follows:
q
r(q)
q
The maximally unstable mode
q
is given by the solution to
r
0
(
q
) = 0, which
we can easily check to be given by
q
=
r
f
00
2κ
.
Now consider the equal time non-equilibrium structure factor
S
q
(t) = hφ
q
(t)φ
q
(t)i S
q
(0)e
2r(q)t
.
As time evolves, this gets more and more peaked around q = q
:
q
S
q
(t)
q
So we see a growth of random structure with scale
L π/q
. This process is
called spinodal decomposition.
L
Note that this computation was done on the assumption that
˜
φ
is small, where
we ignored the quadratic terms. At intermediate
t
, as these phase separated
states grow, the quartic terms are going to kick in. An informal use of variational
theory says we should replace f
00
by
¯
f
00
, where
¯
f
00
= f
00
+
3b
(2π)
4
Z
q
max
S
q
(t) d
d
q.
This says
¯
f
00
is less negative as the fluctuations grow. Since
q
=
r
¯
f
00
2k
,
this moves to a smaller
q
. So
L
(
t
)
π/q
(
t
) starts increasing. This is called
domain growth.
In the late stages, we have large regions of
φ ±φ
B
, so it is almost in
equilibrium locally. We are well away from the exponential growth regime, and
the driving force for domain growth is the reduction of interfacial area. We can
estimate the free energy (per unit volume) as
F
V
=
σA(t)
V
,
where
A
(
t
) is the area of the interface. So by dimensional analysis, this is
σ/L(t). We have calculated the interfacial surface tension σ before to be
σ =
8κa
3
9b
2
1/2
,
but it doesn’t really matter what this is. The ultimate configuration with
minimal surface area is when we have complete phase separation. The result is
that
L(t)
M
σ
φ
B
t
1/3
.
We will not derive this yet, because this result is shared with the late time
behaviour of the third regime, and we will discuss this at that point.
Regime 3
Finally, consider the third regime. Suppose we have
¯
φ
=
φ
B
+
δ
, where
δ
is
small. The system is locally stable, so
r
(
q
)
>
0 for all
q
. On the other hand, it is
globally unstable, so phase separation is preferred. To achieve phase separation,
we must overcome a nucleation barrier, and we must rely on noise to do that.
To understand how the process will look like, formally, we can inspect the
path probabilities
P[φ(r, t)] = N exp
β
4M
Z
|J + Mµ|
2
dr dt
given by the Langevin equation. We seek to find the most likely trajectory from
the initial to the final state. In field theory, this is called the instanton path,
and in statistical physics, this is called large deviation theory. Instead of doing
this, we use our physical intuition to guide ourselves.
Heuristically, we expect that if we start with a uniform phase
φ
=
φ
B
+
δ
,
then at some point, there will be some random small droplet with
φ
= +
φ
B
of
small radius
R
. This is already unlikely, but after this, we need
R
to increase
until we have full phase separation. The key question is how unlikely is this
process?
The idea is to consider the cost of having a droplet of radius
R
. First there is
the cost of having an interface, which is 4
πσR
2
. However, the fact that we have
+
φ
B
areas is energetically favorable, and this grows as the volume. So we get a
cubic term R
3
. If we add these two together, we get a barrier:
R
F (R)
R
F
Once
R > R
, it is then energetically favorable for the radius to continue
increasing, and then we can easily reach phase separation. To reach this, we
must rely on noise to push us over this barrier, and this noise-induced rate is
e
βF
. To see what happens afterwards, we need to better understand how
droplets work.
Droplet in equilibrium
The mechanics of a droplet is slightly less straightforward than what one might
hope, due to the presence of surface tension that tries to compress the droplet.
The result is that the value of
φ
inside and outside the droplet is not exactly
±φ
B
, but with a shift.
For simplicity, we shall first consider an equilibrium system with a droplet.
This is achieved by having a large box of fluid with
¯
φ
just slightly above
φ
B
.
Then in the phase separated state, the +
φ
B
phase will lump together in a droplet
(if we had
¯
φ = 0, then we would have a horizontal interface).
2
1
Within each region 1 and 2, the value of
φ
is constant, so the term that contributes
to the free energy is
f(φ) =
a
2
φ
2
+
b
4
φ
4
.
We would expect 1 and 2 to respectively be located at
φ
f
1 2
When we have a spherical interface, 1 and 2 are not exactly at
±φ
B
. To see this,
Consider the bulk chemical potential
µ =
f
φ
,
The thermodynamic pressure is then
Π = µφ f.
This is the negative of the y-intercept of the tangent line at φ.
If we have a flat interface, which we can think of as the limit
R
, then
we require
µ
bulk
1
= µ
bulk
2
, Π
bulk
1
= Π
bulk
2
.
This means the points 1, 2 have a common tangent
φ
f
1 2
If we have a droplet, then there is surface tension. Consider an imaginary
interface between the upper and lower interface. Then the pressure difference
tries to move the upper hemisphere up, and contributes a force of
2
Π
1
)
πR
2
,
while the interfacial tension pulls the boundary down by 2
π
. In general, in
d
dimensions, we require
Π
2
= Π
1
+
σ
R
(d 1)
This is called the Laplace pressure.
In static equilibrium, we still require
µ
1
=
µ
2
, since this is the same as saying
J
=
µ
= 0. So
f
has the same slope at
φ
1
and
φ
2
. However, the two tangent
lines no longer have a common intercept, but they are separated by
σ
R
(
d
1).
So it looks like
φ
f
φ
1
φ
2
Π
2
Π
1
To solve for this, we take the approximation that
δ
is small for
R
decently
large. Then we can write
f
1
= f(φ
B
+ δ
1
)
1
2
f
00
(φ
B
)δ
2
1
+ f (φ
B
)
f
2
= f(+φ
B
+ δ
2
)
1
2
f
00
(+φ
B
)δ
2
1
+ f (+φ
B
).
So µ
i
= αδ
i
, where α = f
00
(±φ
B
). So we find that up to first order, δ
1
= δ
2
.
To compute δ, we compute
Π
1
= µ
1
φ
1
f
1
= αδφ
B
.
Similarly, we have Π
2
= +αδφ
B
. So
Π
1
Π
2
= 2αφ
B
δ.
Since this equals (d 1)
σ
R
, we have
δ =
d 1
2αφ
B
·
σ
R
.
Multiple droplet dynamics
We now move on to understand multiple droplet dynamics. This is relevant
because we expect that noise will produce multiple droplets around the fluid,
which will then interact and combine to a single phase separated state.
The way droplets interact with each other is that once we have a droplet
of large
φ
, then the average
φ
outside of the droplet will decrease. So to begin
understanding this scenario, we first see how a droplet reacts when the relative
density of the outside bath is not what it expects.
So suppose we have a (3D) droplet of radius
R
inside a bath with
φ
=
φ
B
+
ε
,
where
ε 6
=
δ
=
δ
(
R
). This
ε
is called the supersaturation. Note that to have a
droplet of radius
R
, the value of
φ
inside and immediately outside the droplet
must be
±φ
B
+
δ
. Outside of the droplet, the value of
φ
will slowly decay to
φ
B
+ ε. Thus, outside of the droplet, we write
φ(r) = φ
B
+
˜
φ(r),
where
˜
φ() = ε and
˜
φ(R
+
) = δ.
In this situation, unless
δ
happens to be
ε
, we have a gradient of
φ
, hence a
gradient of chemical potential, hence a flux. Again in Model B, we have
˙
φ = −∇ · J, J = Mµ = M α
˜
φ(r),
assuming a weak enough gradient. We assume this has a quasi-static behaviour,
which is reasonable since molecules move quickly relative to how quickly the
droplet changes in size. So to solve for
˜
φ
(
r
) at any point in time, we set
˙
φ
= 0.
So
2
˜
φ = 0. We solve this with boundary conditions
˜
φ() = ε,
˜
φ(R
+
) = δ.
So we have
˜
φ = ε + (δ ε)
R
r
.
Now if we assume this is what
˜
φ
(
r
) looks like, then the current just outside the
droplet gives
J(R
+
) = Mµ = αM
˜
φ
r
R
+
= αM(δ ε)
B
r
2
r=R
+
=
αM(δ ε)
R
.
Thus, when
δ
and
ε
are not the same, there is a flow of fluid in or out of the
droplet. The discontinuity in
φ
across the boundary is
φ
= 2
φ
B
. So mass
conservation implies
2φ
B
˙
R = J =
αM(δ ε)
R
.
Thus, we conclude that
˙
R =
1
2φ
B
αM
R
(ε δ(R))
.
We can plug in our previous expression for
δ
. Fixing
ε
, we can plot
˙
R
as follows:
R
˙
R
R
where
R
=
σ
αεφ
B
.
So if we have a bath containing many droplets, then the big droplets grow and
the small droplets shrink. Indeed, the interfacial tension penalizes small droplets
more heavily than large droplets.
To understand exactly how these grow, we make a scaling ansatz that there
is only one length scale, namely the mean droplet size
¯
R. Then we have
˙
¯
R
1
2φ
B
αM
¯
R
(ε δ(
¯
R)).
We know that the value of
ε
is also determined by
¯
R
, so we know
ε δ
(
¯
R
) is of
order δ(
¯
R). Hence
˙
¯
R
Mσ
φ
2
B
¯
R
2
So
¯
R
3
Mσt
φ
2
B
.
So the typical droplet size is t
1/3
. Likewise, R
t
1/3
, and so ε t
1/3
.
So if we have a sea of droplets, they go into this competitive process, and
we get fewer and fewer droplets of larger and larger size. This is called Ostwald
ripening, and is a diffusive coarsening mechanism.
We have the same scaling for non-droplet geometries, e.g. spinodal decom-
position at late times. In this case, our domains flatten and enlarge, and we
have
L(t)
Mσ
φ
2
B
t
1/3
.
In practice, we often want to stop this from happening. One way to do so is
to add trapped species insoluble in the continuous phase, e.g. polymers or salt.
If there are
N
particles inside the droplet exerting ideal gas pressure, then we
have
Π
2
Π
1
=
2σ
R
+
Nk
B
T
4/3πR
3
,
We again have µ
1
= µ
2
. This ends up giving a new boundary condition at R
+
,
˜
φ(R
+
) =
σ
αRφ
B
3Nk
B
T
8αφ
B
πR
3
=
1
2αφ
B
Lap
Π
sol
)
The first term is the Laplace pressure just as before, while the second term is
the extra term from the trapped species.
If we put this back into the system of equations before, then we have a new
equation of motion
˙
R =
1
2φ
B
αM
R
ε
σ
αφ
B
R
+
3Nk
B
T
8αφ
B
πR
3

.
R
˙
R
R
s
We now see that there is a stable fixed point
R
s
, and further shrinkage is
prevented by the presence of the trapped species that will be further compressed
by shrinking. Thus, if we manage to produce a system with all droplets of size
< R
, then we end up with a lot of small but finite size droplets R
s
.
7 Model H
Model B was purely diffusive, and the only way
φ
can change is by diffusion.
Often, in real life, fluid can flow as well. If the fluid has a velocity
v
, then our
equation is now
˙
φ + v · φ = −∇ · J.
The
v · φ
is called the advection term. Our current
J
is the same as before,
with
J = M
δF
δφ
+
p
2k
B
T MΛ.
We also need an evolution equation for
v
, which will be the Navier–Stokes
equation with some noise term. We assume flow is incompressible, so
· v
= 0.
We then have the Cauchy equation with stress tensor Σ
TOT
, given by
ρ(
˙
v + v · v) = · Σ
TOT
+ body forces.
We will assume there is no body force. This is essentially the momentum
conservation equation, where Σ
TOT
is the momentum flux tensor.
Of course, this description is useless if we don’t know what Σ
TOT
looks like.
It is a sum of four contributions:
Σ
TOT
= Σ
p
+ Σ
η
+ Σ
φ
+ Σ
N
.
The Σ
p
term is the pressure term, given by
Σ
p
ij
= P δ
ij
.
We should think of this P as a Lagrange multiplier for incompressibility.
The Σ
η
term is the viscous stress, which we can write as
Σ
η
ij
= η(
i
v
j
+
j
v
i
)
For simplicity, we assume we have a constant viscosity
η
. In general, it
could be a function of the composition.
The Σ
φ
term is the φ-stress, given by
Σ
φ
= Πδ
ij
κ(
i
φ)(
j
φ), Π = φµ F.
This is engineered so that
· Σ
φ
= φµ.
This says a non-constant chemical potential causes things to move to even
that out.
The final term is a noise stress with
hΣ
N
ij
(r, t
N
k`
(r
0
, t
0
)i = 2k
B
T η
δ
ik
δ
j`
+ δ
i`
δ
jk
2
3
δ
ij
δ
k`
δ(rr
0
)δ(tt
0
).
The last term
δ
ij
δ
k`
is there to ensure the noise does not cause any
compression. This is a white noise term whose variance is determined by
the fluctuation dissipation theorem.
We can then compute
· Σ
TOT
= · Σ
P
+ · Σ
η
+ · Σ
φ
+ · Σ
N
= −∇P + η
2
v + φµ + · Σ
N
Hence Model H has equations
˙
φ + v · φ = −∇ · J
J = M µ +
p
2k
B
T MΛ
· v = 0
ρ(
˙
v + v · v) = η
2
v P φµ + · Σ
N
.
Compared to Model B, we have the following new features:
(i) φµ drives deterministic fluid flow.
(ii) Σ
N
drives a random flow.
(iii) Fluid flow advects φ.
How does this affect the coarsening dynamics? We will see that (i) and (iii) gives
us enhanced coarsening of bicontinuous states. However, this does not have any
effect on isolated/disconnected droplet states, since in a spherically symmetric
setting,
φµ
and
P
will be radial, and so
· v
= 0 implies
v
= 0. In other
words, for φµ to drive a flow, we must have some symmetry breaking.
Of course, this symmetry breaking is provided by the noise term in (ii). The
result is that the droplets will undergo Brownian motion with
hr
2
i Dt
, where
D =
k
B
T
4πηR
is the diffusion constant.
If we think about the Ostwald process, even if we manage to stop the small
droplets from shrinking, they may collide and combine to form larger droplets.
This forms a new channel for instability, and separate measures are needed to
prevent this. For example, we can put charged surfactants that prevent collisions.
We can roughly estimate the time scale of this process. We again assume
there is one length scale
¯
R
(
t
), which determines the size and separation of
droplets. We can then calculate the collision time
t '
¯
R
2
D(
¯
R)
¯
R
3
η
k
B
T
.
Each collision doubles the volume, and so
¯
R
2
1/3
¯
R
. Taking the logarithm, we
have
log
¯
R
log 2
3
in time t.
So we crudely have
log
¯
R
t
log 2
3
k
B
T
η
¯
R
3
.
If we read this as a differential equation, then we get
d log
¯
R
dt
=
1
¯
R
˙
¯
R
k
B
T
η
¯
R
3
.
So we find that
¯
R
2
˙
¯
R
k
B
T
η
.
So
¯
R(t)
k
B
T
η
t
1/3
.
This is diffusion limited coalescence.
Recall that in the Ostwald process, droplets grew by diffusion of molecules,
and we had the same power of t. However, the coefficient was different, with
¯
R
Mσ
φ
B
t
1/3
.
It makes sense that they have the same scaling law, because ultimately, we are
still doing diffusion on different scales.
Bicontinuous states
We now see what the fluid flow does to bicontinuous states.
L
Again assume we have a single length scale
L
(
t
), given by the domain size.
Again we assume we have a single length scale. As time goes on, we expect
L
(
t
)
to increase with time.
A significant factor in the evolution of the bicontinuous phase is the Laplace
pressure, which is ultimately due to the curvature
K
. Since there is only one
length scale, we must have
˙
K v.
The Laplace pressure then scales as
σ
L
.
The noise terms Σ
N
matter at early times only. At late times, the domains
grow deterministically from random initial conditions. The key question is how
L(t) scales with time. The equation of motion of the flow is
ρ(
˙
v + v · v) = η
2
v P φµ,
We make single length scale approximations as before, so that
v
˙
L
and
L
1
. Then we have an equation of the form.
ρ
¨
L + ρ
˙
L
2
L
η
˙
K
L
2
+ Lagrange multiplier +
σ
L
2
()
where we recall that at curved interfaces,
µ ±
σ
R
. Here we have a single variable
L
(
t
), and three dimensionful parameters
ρ, η, σ
. We can do some dimensional
analysis. In d dimensions, we have
ρ = M L
d
, η = ML
2d
T
1
, σ = ML
3d
T
2
.
We want to come up with combinations of these for length to depend on time,
and we find that in three dimensions, we have
L
0
=
η
2
ρσ
, t
0
=
η
3
ρσ
2
.
One can check that these are the only combinations with units
L, T
. So we must
have
L(t)
L
0
= f
t
t
0
.
We now substitute this into the equation (
). We then get a non-dimensionalized
equation for ()
αf
00
+ βf
02
/f = γ
f
0
f
2
+
δ
f
2
,
with α, β, γ, δ = O(1) dimensionless numbers.
If we think about this, we see there are two regimes,
(i) The LHS (inertia) is negligible at small t/t
0
(or small f ). Then we get
γf
0
f
2
+
δ
f
2
= 0,
so
f
0
is a constant, and so
L
grows linearly with
t
. Putting all the
appropriate constants in, we get
L
σ
η
t.
This is called the viscous hydrodynamic regime, VH .
(ii)
For large
f
, we assume we have a power law
f
(
x
)
x
y
, where
y >
0 (or
else f would not be large). Then
¯αx
y2
+
¯
βx
y2
= ¯γx
y1
+ δx
2y
.
It turns out at large
f
, the
x
y1
term is negligible, scaling wise. So we
have y 2 = 2y, or equivalently, y =
2
3
. So
L
L
0
t
t
0
2/3
.
Putting back our factors, we have
L
σ
ρ
1/3
t
2/3
.
This is called the inertial hydrodynamic regime, IH . In this regime, interfa-
cial energy is converted into kinetic energy, and then only “later” dissipated
by η.
Essentially, in the first regime, the system is overdamped. This happens until
the right-hand side becomes big and takes over, until the viscous term finally
takes over. In practice, it is difficult to reach the last regime in a lab, since the
time it takes is often 10
4
.
Droplet vs bicontinuous
In general, when do we expect a bicontinuous phase and when do we expect a
droplet phase?
In three dimensions, the rule of thumb is that if
ψ
0
.
4
to
0
.
6, then
we always get a bicontinuous medium. If
ψ <
0
.
3 or
>
0
.
7, then we
get droplets always. In between these two regions, we initially have
bicontinuous medium, which essentially de-percolates into droplets.
In two dimensions, things are different. In the fully symmetric case, with a
constant
η
throughout, and
F
is strictly symmetric (
F
(
φ
) =
F
(
φ
)), the
only case where we have bicontinuous phase is on the ψ =
1
2
line.
8 Liquid crystals hydrodynamics
8.1 Liquid crystal models
We finally turn to the case of liquid crystals. In this case, our order parameter no
longer takes a scalar value, and interesting topological phenomena can happen.
We first write down the theory in detail, and then discuss coarsening behaviour
of liquid crystals. We will describe the coarsening purely via geometric means,
and the details of the model are not exactly that important.
Recall that there are two types of liquid crystals:
(i)
Polar liquid crystals, where the molecules have “direction”, and the order
parameter is p, a vector, which is orientational but not positional.
(ii)
Nematic liquid crystals, where the molecules are like a cylinder, and there
is an order parameter Θ.
We will do the polar case only, and just quote the answers for the nematic
case.
As before, we have to start with a free energy
F =
Z
a
2
|p|
2
+
b
4
|p|
4
+
κ
2
(
i
p
j
)(
i
p
j
)
dr
Z
F dr.
The first two terms are needed for the isotropic-polar transition to occur. Note
that
(i) F [p] = F [p], so we have no cubic term.
(ii)
A linear term would correspond to the presence of an external field, e.g.
magnetic field.
(iii)
The
κ
term penalizes splay, twist and bend, and this term penalizes them
roughly equally. This is called the one elastic constant approximation. If
we want to be more general, we need something like
κ
1
2
|∇ · p|
2
+
κ
2
2
|
ˆ
p · p|
2
+
κ
3
2
|
ˆ
p p|
2
.
Here
p
is not conserved, so
˙
p
is not of the form
−∇ · J
. Instead, (without flow)
we have
˙
p = Γh, h =
δF
δp(r)
,
where
h
(
r
) is called the molecular field, and Γ is a constant, called the angular
mobility.
We now want to generalize this to the case when there is a field flow. We
can just write this as
Dp
Dt
= Γh,
where D is some sort of comoving derivative. For the scalar field, we had
Dφ
Dt
=
˙
φ + v · φ.
Note that the advective term is trilinear, being first order in
v
,
and
φ
. For
p
,
there is certainly something like this going on. If we have a translation, then
p
gets incremented by
p = v · p t,
as for a scalar.
There is at least one more thing we should think about, namely if
v
is
rotational, then we would expect this to rotate
p
as well. We have the corotational
term
p = ω p t,
where ω is the angular velocity of the fluid, given by
ω
i
=
1
2
ε
ijk
jk
,
jk
=
1
2
(
i
v
j
j
v
i
).
This part must be present. Indeed, if we rotate the whole system as a rigid body,
with v(r) = ω × r, then we must have
˙
p = ω p.
It turns out in general, there is one more contribution to the advection, given
by
p = ξD · pt
with D
ij
=
1
2
(
i
v
j
+
j
v
i
) and
ξ
a parameter. This is the irrotational part of
the derivative. The reason is that in a general flow,
p
needn’t simply rotate with
ω. Instead, it typically aligns along streamlines. In total, we have
Dp
Dt
= (
t
+ v · )p + · p ξD · p.
The parameter
ξ
is a molecular parameter, which depends on the liquid crystal.
The ξ = 1 case is called no slip, and the ξ = 0 case is called full slip.
With this understanding, we can write the hydrodynamic equation as
Dp
Dt
= Γh, h =
δF
δp
.
We next need an equation of motion for v. We can simply
ρ(
t
+ v · )v = η
2
v P + · Σ
p
,
where Σ
p
is a stress tensor coming from the order parameter.
To figure out what Σ
p
should be, we can consider what happens when we
have an “advective” elastic distortion. In this case, we have
Dp
Dt
= 0, so we have
˙
p = v · p · p + ξD · p,
The free energy change is then
δF =
Z
δF
δp
·
˙
p t dr = ∆t
Z
h · p dr,
On the other hand, the free energy change must also be given by
δF =
Z
Σ
p
ij
i
u
j
(r) dr,
the product of the stress and strain tensors. By matching these two expressions,
we see that we can write
Σ
p
ij
= Σ
(1)
ij
+ Σ
(2)
ij
+ Σ
(3)
ij
,
where
i
Σ
(1)
ij
= p
k
j
h
k
, Σ
(2)
ij
=
1
2
(p
i
h
j
p
j
h
i
), Σ
(2)
ij
=
ξ
2
(p
i
h
j
+ p
j
h
i
).
Analogous results for a nematic liquid crystal holds, whose derivation is a
significant pain. We have
DQ
Dt
= ΓH, H
ij
=
δF
δQ
ij
Tr
δF
δQ
δ
ij
d
.
The second term in H is required to ensure Q remains traceless.
The total derivative is given by
DQ
Dt
= (
t
+v ·)Q+ (Ω ·Q Q ·ω) + ξ(D ·Q + Q ·D)2ξ
Q +
1
d
Tr(Q ·v).
The terms are the usual advection term, rotation, alignment/slip and tracelessness
terms respectively. The Navier–Stokes equation involves a stress term
Σ
Q
= Σ
Q,1
+ Σ
Q,2
,
where
k
Σ
Q,1
k,`
= Q
ij
`
H
ij
Σ
Q,2
= Q · H H · Q ξ
2
3
H + 2
d
QH 2Q tr(QH)
,
with the hat denoting the traceless symmetric part. The rest of the Navier–Stokes
equations is the same.
8.2 Coarsening dynamics for nematics
We will discuss the coarsening dynamics for nematic liquid crystals, and indicate
how polar liquid crystals are different when appropriate. As before, we begin
with a completely disordered phase with
Q
= 0, and then quench the system
by dropping the temperature quickly. The liquid crystals then want to arrange
themselves. Afterwards, we local have
Q =
λ 0 0
0 λ/2 0
0 0 λ/2
with free energy
f
(
λ
). If the quench is deep enough, we have spinodal-like
instability, and we quickly get locally ordered. Coordinate-independently, we
can write
Q =
ˆ
λ
n
i
n
j
1
3
δ
ij
.
Since all this ordering is done locally, the principal axis
n
can vary over space.
There is then a slower process that sorts out global ordering, driven by the elastic
part of the free energy.
Compare this with Model B/H: at early times, we have
φ
:
±φ
B
, but we get
domain walls between the ±φ
B
phases.
x
φ
The late time dynamics is governed by the slow coarsening of these domain walls.
The key property of this is that it has reduced dimensionality, i.e. the domain
wall is 2 dimensional while space is 3 dimensional, and it connects two different
grounds states (i.e. minima of F ).
The domain wall can be moved around, but there is no local change that
allows us to remove it. They can only be removed by “collision” of domain walls.
Analogous structures are present for nematic liquid crystals. The discussion
will largely involve us drawing pictures. For now, we will not do this completely
rigorously, and simply rely on our intuitive understanding of when defects can
or cannot arise. We later formalize these notions in terms of homotopy groups.
We first do this in two dimensions, where defects have dimension
<
2. There
can be no line defects like a domain wall. The reason is that if we try to construct
a domain wall
then this can relax locally to become
On the other hand, we can have point defects, which are 0-dimensional. Two
basic ones are as follows:
q =
1
2
q = +
1
2
The charge
q
can be described as follows we start at a point near the
defect, and go around the defect once. When doing so, the direction of the order
parameter turns. After going around the defect once, in the
q
=
1
2
, the order
parameter made a half turn in the opposite sense to how we moved around the
defect. In the q = +
1
2
case, they turned in the same sense.
We see that
q ±
1
2
are the smallest possible topological charge, and is a
quantum of a charge. In general, we can have defects of other charges. For
example, here are two q = +1 charges:
hedgehog
vortex
Both of these are
q
= +1 defects, and they can be continuously deformed
into each other, simply by rotating each bar by 90
. For polar liquid crystals,
the quantum of a charge is 1.
If we have defects of charge greater than
±
1
2
, then they tend to dissociate
into multiple
q
=
±
1
2
defects. This is due to energetic reasons. The elastic
energy is given by
F
ell
=
κ
2
|∇ · Q|
2
κ
2
λ
2
|( · n)n + n · n|
2
.
If we double the charge, we double the
Q
tensor. Since this term is quadratic
in the gradient, putting two defects together doubles the energy. In general,
topological defects tend to dissociate to smaller q-values.
To recap, after quenching, at early stages, we locally have
Q 2λ(nn
1
2
1).
This
n
(
r
) is random, and tend to vary continuously. However, topological defects
are present, which cannot be ironed out locally. All topological defects with
|q| >
1
2
dissociate quickly, and we are left with q = ±
1
2
defects floating around.
We then have a late stage process where opposite charges attract and then
annihilate. So the system becomes more and more ordered as a nematic. We
can estimate the energy of an isolated defect as
˜κ
2
|( · n)n + n|
2
,
where ˜κ = κλ
2
. Dimensionally, we have
1
r
.
So we have an energy
E ˜κ
Z
1
r
2
dr ' ˜κ log
L
r
0
,
where
L
is the mean spacing and
r
0
is some kind of core radius of the defect.
The core radius reflects the fact that as we zoom close enough to the core of the
signularity,
λ
is no longer constant and our above energy estimate fails. In fact,
λ 0 at the core.
Recall that the electrostatic energy in two dimensions is given by a similar
equation. Thus, this energy is Coulombic, with force
˜κ
R
. Under this force, the
defects move with overdamped motion, with the velocity being proportional to
the force. So
˙
R
1
R
,
˙
L
1
L
.
So
L(t) t
1/2
.
This is the scaling law for nematic defect coarsening in 2 dimensions.
8.3 Topological defects in three dimensions
In three dimensions, we also have defects of the above kind lying along a line.
For such line defects, everything we said so far carries through separated at a
distance R, the interaction force is
˜κ
R
and so and so we similarly have
L(t) t
1/2
.
However, in three dimensions, the
q
=
±
1
2
defects are the same topologically. In
other words, we can change +
q
to
q
via continuous, local deformations. This
involves rotating out to the
z
direction, which is not available in two dimensions.
While it is possible to visually understand how this works, it is difficult to draw
on paper, and it is also evident that we should proceed in more formal manners
to ensure we understand exactly how these things work.
To begin, we have a space
M
of order parameters. In our case, this is the
space of all possible orientations of rods.
Example.
In the case of a polar liquid crystal in
d
dimensions, we have
M
=
S
d1
, the (d 1)-dimensional unit sphere.
Example.
For nematic liquid crystals in
d
-dimensions, we have
M
=
RP
d1
,
which is obtained from S
d1
by identifying antipodal points.
When we discussed the charge of a topological defect, we picked a loop around
the singularity and see what happened when we went around a defect. So we
pick a domain
D
that encloses a defect core, and consider the map
f
:
D M
that assigns to each point the order parameter at that point. In our cases,
D
is
a circle S
1
, and so f is a loop in M.
We say two mappings
f
1
, f
2
, are homotopic if they can be continuously
deformed into each other. Defects lie in the same homotopy class if maps for
all
D
’s enclosing them are homotopic. The fundamental group
π
1
(
M
) is the set
of all homotopy classes of maps
S
1
M
. This encodes the set of all possible
charges.
Since we call it a fundamental group, it had better have a group structures.
If we have two defects, we can put them next to each other, and pick a new
circle that goes around the two defects. This then gives rise to a new homotopy
class S
1
M.
More generally, if we consider
d n
-dimensional defects, then we can enclose
the defect with a sphere of dimension
n
1. The corresponding classes live in
the higher homotopy groups π
n1
(M).
Example.
Observe that
RP
1
is actually just
S
1
in disguise, and so
π
1
(
RP
1
) =
Z
.
The generator of π
1
(RP
1
) is the charge
1
2
topological defect.
Example. We can visualize RP
2
as a certain quotient of the disk, namely
where we identify the two arcs in the boundary according to the arrow. Observe
that the two marked points are in fact the same point under the identification.
If we have a path from the first point to the second point, then this would be
considered a loop in RP
2
, and this is the q =
1
2
defect.
Observe that in the two-dimensional case, the
q
=
±
1
2
defects correspond to
going along the top arc and bottom arc from the left to right respectively. In
RP
2
, there is then a homotopy between these two paths by going through the
disk. So in RP
2
, they lie in the same homotopy class.
In general, it is easy to see that
π
1
(
RP
2
) =
Z/
2
Z
, so
q
=
1
2
is the unique
non-trivial defect.
This is particularly interesting, because two
q
=
1
2
defects can merge and
disappear! Similarly, what you would expect to be a
q
= 1 defect could locally
relax to become trivial.
Observe that in our “line defects”, the core can actually form a loop instead.
We can also have point defects that correspond to elements in
π
2
(
M
)
=
Z
. It is
an exercise to draw some pictures yourself to see how these look.
9 Active Soft Matter
We shall finish the course by thinking a bit about motile particles. These are
particles that are self-propelled. For example, micro-organisms such as bacteria
and algae can move by themselves. There are also synthetic microswimmers. For
example, we can make a sphere and coat it with gold and platinum on two sides
PtAu
We put this in hydrogen peroxide H
2
O
2
. Platinum is a catalyst of the decompo-
sition
2H
2
O
2
2H
2
O + O
2
,
and this reaction will cause the swimmer to propel forward in a certain direction.
This reaction implies that entropy is constantly being produced, and this cannot
be captured adequately by Newtonian or Lagrangian mechanics on a macroscopic
scale.
Two key consequences of this are:
(i) There is no Boltzmann distribution.
(ii)
The principle of detailed balance, which is a manifestation of time reversal
symmetry, no longer holds.
Example. Take bacteria in microfluidic enclosure with funnel gates:
In this case, we expect there to be a rotation of particles if they are self-
propelled, since it is easier to get through one direction than the other. Contrast
this with the fact that there is no current in the steady state for any thermal
equilibrium system. The difference is that Brownian motion has independent
increments, but self-propelled particles tend to keep moving in the same direction.
Note also that we have to have to break spatial symmetry for this to happen.
This is an example of the Ratchet theorem, namely if we have broken time
reversal symmetry pathwise, and broken spatial symmetry, then we can have
non-zero current.
If we want to address this type of system in the language we have been using,
we need to rebuild our model of statistical physics. In general, there are two
model building strategies:
(i)
Explicit coarse-graining of “micro” model, where we coarse-grain particles
and rules to PDEs for ρ, φ, P, Q.
(ii)
Start with models of passive soft matter (e.g. Model B and Model H), and
add minimal terms to explicitly break time reversal phenomenologically.
Of course, we are going to proceed phenomenologically.
Active Model B
Start with Model B, which has a diffusive, symmetric scalar field
φ
with phase
separation:
˙
φ = −∇ · J
J = −∇˜µ +
2DΛ.
We took
F =
Z
a
2
+
b
4
φ
4
+
κ
2
(φ)
2
dr.
To model our system without time reversal symmetry, we put
˜µ =
δF
δφ
+ λ(φ)
2
.
The new term breaks the time reversal structure. These equations are called
active Model B. Another way to destroy time reversal symmetry is by replacing
the white noise with something else, but that is complicated
Note that
(i)
(
φ
)
2
is not the functional derivative of any
F
. This breaks the free energy
structure, and
P
F
P
B
6= e
β(F
2
F
1
)
for any F [φ]. So time reversal symmetric is broken barring miracles.
(ii)
We cannot achieve the breaking by introducing a polynomial term, since if
g(φ) is a polynomial, then
g(φ) =
δ
δφ
Z
dr
Z
φ
g(u) du
!
.
So gradient terms are required to break time reversal symmetry. We will
later see this is not the case for liquid crystals.
(iii)
The active model B is agnostic about the cause of phase separation at
a < 0. There are two possibilities:
(a) We can have attractive interactions
(b)
We can have repulsive interactions plus motile particles: if two parti-
cles collide head-on, then we have pairwise jamming. They then move
together for a while, and this impersonates attraction. This is called
MIPS mobility-induced phase separation. It is possible to study
this at a particle level, but we shall not.
(iv)
The dynamics of coarsening during phase separation turns out to be similar,
with L(t) t
1/3
. The Ostwald–like process remains intact.
(v)
The coexistence conditions are altered. We previously found the coexistence
conditions simply by global free energy minimization. In the active case,
we can’t do free energy minimization, but we can still solve the equations
of motion explicitly. In this case, instead of requiring a common tangent,
we have equal slope but different intercepts, where we set
(µφ f)
1
= (µφ f)
2
+ .
This is found by solving J = 0, so
˜µ =
f
φ
κ
2
φ + λ(φ)
2
= const.
(vi) There is a further extension, active model B+, where we put
J = −∇˜µ +
2DΛ + ζ(
2
φ)φ.
This extra term is similar to
(
λ
(
φ
)
2
) in that it has two
φ
’s and three
’s, and they are equivalent in 1 dimension, but in 1 dimension only. This
changes the coarsening process significantly. For example, Ostwald can
stop at finite R (see arXiv:1801.07687).
Active polar liquid crystals
Consider first a polar system. Then the order parameter is
p
. In the simplest
case, the field is relaxational with v = 0. The hydrodynamic level equation is
˙
p = Γh, h =
δF
δp
.
We had a free energy
F =
Z
a
2
p
2
+
b
4
p
4
+
κ
2
(
α
p
β
)(
α
p
β
)
dr.
As for active model B,
h
can acquire gradient terms that are incompatible with
F
.
But also, we can have a lower order term in
that is physically well-motivated —
if we think of our rod as having a direction
p
, then it is natural that
p
wants to
translate along its own direction at some speed
w
. Thus,
p
acquires self-advected
motion wp. Thus, our equation of motion becomes
˙
p + p · p = Γh.
This is a bit like the Navier–Stokes equation non-linearity. Now couple this to a
fluid flow v. Then
Dp
Dt
= Γh,
where
Dp
Dt
=
t
+ v ·
p + · p ξD · p + wp · p.
The Navier–Stokes/Cauchy equation is now
(
t
+ v · )v = η
2
v P + · Σ
(p)
+ · Σ
A
,
where as before,
· Σ
(p)
= p
i
j
h
j
+
i
1
2
(p
i
h
j
p
j
h
i
) +
ξ
2
(p
i
h
j
+ p
j
h
i
)
.
and we have a new term Σ
A
given by the active stress, and the lowest order
term is
ζp
i
p
j
. This is a new mechanical term that is incompatible with
F
. We
then have
· Σ
A
= ( · p)p.
We can think of this as an effective body force in the Navier–Stokes equation.
The effect is that we have forces whenever we have situations looking like
or
In these cases, We have a force acting to the right for
ζ >
0, and to the left if
ζ < 0.
These new terms give spontaneous flow, symmetric breaking and macroscopic
fluxes. At high w, ζ, we get chaos and turbulence.
Active nematic liquid crystals
In the nematic case, there is no self-advection. So we can’t make a velocity from
Q. We again have
DQ
Dt
= ΓH, H =
δF
δQ
traceless
.
where
DQ
Dt
is given by
DQ
Dt
= (
t
+ v · )Q + S(Q, K, ξ).
Here K = v and
S = ( · Q Q · Ω) ξ(D · Q + Q · D) + 2ξ
Q +
1
d
Tr(Q · K)
Since there is no directionality as in the previous case, the material derivative
will remain unchanged with active matter. Thus, at lowest order, all the self-
propelled motion can do is to introduce an active stress term. The leading-order
stress is
Σ
A
= ζQ.
This breaks the free energy structure. Indeed, if we have a uniform nematic, then
the passive stress vanishes, because there is no elastic distortion at all. However,
the active stress does not since
ζQ 6
= 0. Physically, the non-zero stress is due to
the fact that the rods tend to generate local flow fields around themselves to
propel motion, and these remain even in a uniform phase.
After introducing this, the effective body force density is
f = · Σ
A
= ζ · Q ζλ( · n)n.
This is essentially the same as the case of the polar case. Thus, if we see
something like
then we have a rightward force if ζ > 0 and leftward force if ζ < 0.
This has important physical consequences. If we start with a uniform phase,
then we expect random noise to exist, and then the active stress will destablize
the system. For example, if we start with
and a local deformation happens:
then in the
ζ >
0 case, this will just fall apart. Conversely, bends are destabilized
for
ζ <
0. Either case, there is a tendency for these things to be destabilized,
and a consequence is that active nematics are never stably uniform for large
systems. Typically, we get spontaneous flow.
To understand this more, we can explicitly describe how the activity parameter
ζ affects the local flow patterns. Typically, we have the following two cases:
ζ > 0
contractile
ζ < 0
extensile
Suppose we take an active liquid crystal and put it in a shear flow. A rod-like
object tends to align along the extension axis, at a 45
angle.
If the liquid crystal is active, then we expect the local flows to interact with
the shear flow. Suppose the shear rate is v
x
= yg. Then the viscous stress is
Σ
η
= ηg
0 1
1 0
.
We have
Σ
A
ζλ
nn
1
d
= ζλ
0 1
1 0
if
n
is at 45
exactly. Note that the sign of
ζ
affects whether it reinforces or
weakens the stress.
A crucial property is that Σ
A
does not depend on the shear rate. So in the
contractile case, the total stress looks like
g
Σ
TOT
In the extensile case, however, we have
g
Σ
TOT
g
This is very weird, and leads to spontaneous flow at zero applied stress of
the form
Defect motion in active nematics
For simplicity, work in two dimensions. We have two simple defects as before
q =
1
2
q = +
1
2
Note that the
q
=
1
2
charge is symmetric, and so by symmetry, there cannot be
a net body stress. However, in the
q
= +
1
2
defect, we have a non-zero effective
force density.
So the defects themselves are like quasi-particles that are themselves active.
We see that contractile rods move in the direction of the opening, and the
extensile goes in the other direction. The outcome of this is self-sustaining
“turbulent” motion, with defect pairs
±
1
2
are formed locally. The
1
2
stay put
and the +
1
2
ones self-propel, and depending on exactly how the defect pairs are
formed, the +
1
2
defect will fly away.
Experimental movies of these can be found in T. Sanchez Nature 491, 431
(2012). There are also simulations in T. Shendek, et al, Soft Matter 13, 3853
(2017).