III Theoretical Physics of Soft Condensed Matter (Full)

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

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

∼

–10

Pascals (compare with steel,

which has a shear modulus of 10

Pascals).

–

The bulk modulus

remains large, with order of magnitude

K ∼

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:

−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

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

of quantum fluctuations is given by

~ω

' k

At room temperature, this gives that

∼ 10

−13

, 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

Φ =

πa

where N is the number of colloid particles.

Experimentally, we observe that when Φ

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 =

6πη

where

is the solvent viscosity. Thus, the time

it takes for the particle to

move through a distance of its own radius

is given by

Dτ

, which we can

solve to give

τ ∼

In general, this is much longer than the time scale

of quantum fluctuations,

since a

 ~.

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

In both cases, we see that the elasticity and time scales are given in terms of

. If we ignore thermal fluctuations, then we have

= 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

(

) that

tells us the colloid density near the point

. 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

[

(

)], 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

[

(

)]. 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

. The hydrodynamic PDE

is exactly the Navier–Stokes equation. Assuming incompressibility, so that

˙ρ

= 0,

we get

ρ(

v + v · ∇v) = η∇

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) = η∇

v − ∇p + ∇·Σ

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Σ

(r, t)Σ

, t

)i = 2k

T η(δ

+ δ

)δ(r − r

)δ(t − t

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

log p

where

is Boltzmann’s constant,

is a microstate — a complete specification

of the microscopics (e.g. the list of all particle coordinates and velocities) — and

is the probability of being in a certain microstate.

The axiom of Gibbs is that a system in thermal equilibrium maximizes

subject to applicable constraints.

Example.

In an isolated system, the number of particles

, the energy

and

the volume

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

= 1.

Gibbs says we should maximize

. Writing

for the Lagrange multiplier

maintaining this constraint, we require

∂

∂p

S − λ

= 0.

So we find that

−k

log p

+ 1 − λ = 0

for all i. Thus, we see that all p

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

is no longer fixed, and fluctuates around some average

hEi

. So we can

apply Gibbs’ principle again, where we now sum over all states of all

, with

the restrictions

= 1.

So our equation is

∂

∂p



S − λ

− λ



= 0.

Differentiating this with respect to p

, we get

−k

(log p

+ 1) − λ

− λ

= 0.

So it follows that

−βE

where Z =

−βE

and β = λ

. This is the Boltzmann distribution.

What is this mysterious

? Recall that the Lagrange multiplier

measures

how S reacts to a change in

E. In other words,

∂S

∂E

= λ

= k

β.

Moreover, by definition of temperature, we have

∂S

∂E



V,N,...

So it follows that

β =

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

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

(

)

of space that describes the “average magnetization around

”. 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)] =

i∈ψ

−βE

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[ψ] =

−βF [ψ]

TOT

, Z

TOT

−β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

TOT

= −k

T log Z

TOT

In theory, one can obtain

[

] by explicitly doing a coarse graining of the

macroscopic laws.

Example.

Consider an interacting gas with

particles. We can think of

the energy as a sum of two components, the ideal gas part (

NkT

), and an

interaction part, given by

int

i6j

U(r

− r

where

i, j

range over all particles with positions

, r

respectively, and

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

int

U(r − r

)ρ(r)ρ(r

) dr dr

Similarly, up to a constant, we can write the entropy as

S[ρ] = −k

ρ(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

[

], 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

and

. For

simplicity, we assume we are in the symmetric case, where

and

are the same

“type” of fluids. In other words, the potentials between the fluids are such that

(r) = U

(r) 6= U

(r).

We consider the case where

and

repulse each other (or rather, repulse each

other more than the

and

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

(

) and

(

) be the coarse-grained particle density of each fluid,

and we set our order parameter to be

φ(r) =

(r) − ρ

(r)

+ N

)/V

with

and

, the total amount of fluids

and

, and

the volume. This

is normalized so that φ(r) ∈ [−1, 1].

We model our system with Landau–Ginzburg theory, with free energy given

βF =



| {z }

f(φ)

(∇φ)



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

cφ

term will integrate out to give

, 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

(

∇φ

)

captures at order

∇

(2)

the non-

locality of E

int

i,j∈{A,B}

(r)ρ

(|r − r

|) dr dr

If we assume

(

) is slowly varying on the scale of interactions, then we can

Taylor expand this E

int

and obtain a (∇φ)

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, κ >

while

can take either sign. In general, the sign of

is what determines the

behaviour of the system, so for simplicity, we suppose

and

are fixed, and let

a vary with temperature.

To do mean field theory, we find a single

that minimizes

. Since the

gradient term

(

∇φ

)

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(

φ) =

The global of this function depends only on the sign of

. For

a >

0 and

a <

respectively, the plots look like this:

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 =



f(φ(r)) +

(∇φ)



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

(

) has two minima,

1,2

= ±φ

, where

−a

Now suppose

lies between

±φ

. Then it might be advantageous to have some

parts of the fluid being at

−φ

and the others at

, 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 [−φ

, φ

Suppose there is

many fluid with

, and

many fluid with

Then these quantities must obey

+ V

= V

φ,

+ V

= V.

Concavity tells us we must have

f(φ

) + V

f(φ

) < (V

+ V

)f(

φ).

Thus, if we only consider the

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 → ±∞

, as in our above

scenario. We then find that the thickness of the interface is (of the order)

−2κ

and the cost per unit area of this interface is

σ =



−8κa



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

−1 1

a(T ) = 0

Within the solid lines, we have phase separation, where the ground state of the

system for the given

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

(

)

0, which

we solve to be

−a

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

(

)

to our

βF

? It turns out we can remove

the

(

) term by a linear shift of

and

, 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

, but mod it out by the relation

n ∼ −n

. One

way to do so is to consider the second-rank traceless tensor

. This has the

property that

is the component of a vector

in the direction of

, and

is invariant under

7→ n

. Observe that if we normalize

to be a unit vector,

then

has trace 1. Thus, if we have isotropic rods in

dimensions, then we

have

i =

In general, we can defined a coarse-grained order parameter to be

(r) = hn

local

−

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

is no longer

conserved., i.e. the “total Q”

(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

(

), which has to be a scalar

built on Q. The possible terms are as follows:

(i) There is only one linear one, namely Q

= Tr(Q), but this vanishes.

(ii)

We can construct a quadratic term

(

), and this is in general

non-zero.

(iii)

There is a cubic term

(

), and is also in general non-zero.

(iv) There are two possible quartic terms, namely Tr(Q

)

and Tr(Q

So we can write

f(Q) = a Tr(Q

) + c Tr(Q

) + b

Tr(Q

)

+ b

Tr(Q

This is the local part of the free energy up to fourth order in

. We can go on,

and in certain conditions we have to, but if these coefficients

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

, and are agnostic about the other two

directions, then Q will be given by





−λ/2 0 0

0 −λ/2 0

0 0 λ





, λ > 0.

If the rod is agnostic about the

and

directions, but instead avoids the

-direction, then

takes the same form but with

λ <

0. For the purposes of

(

), we can locally diagonalize

, 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





+ c





+ b





+ b





= ¯aλ

+ ¯cλ

bλ

We can think of this in a way similar to the binary fluid, where

is are sole

order parameter. We fix

and

¯c <

0, and then vary

¯a

. In different situations,

we get

α < α

α = α

α > α

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

∇

= κ

∇

Tr(Q

)

(∇

)(∇

) = κ

(∇ · Q)

(∇

)(∇

) = 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

= 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[φ] =

L(φ, ∇φ) dr.

Under a small change

φ 7→ φ

δφ

(

) with

δφ

= 0 on the boundary, our functional

becomes

A[φ + δφ] =



L(φ, ∇φ) + δφ

∂L

∂φ

+ ∇dφ ·

∂L

∂φ



= A[φ] +

δφ



∂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

by the single variable

, and

by position x. We then have

A =

L(x, ˙x) dt.

Then we have

δA

δx(t)

∂L

∂x

−



∂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 [φ] =



(∇φ)



dr.

Then

δF

δφ(r)

= aφ + bφ

− κ∇

φ.

In mean field theory, we set this to zero, since by definition, we are choosing

a single

(

) that minimizes

. In the first example sheet, we find that the

minimum is given by

φ(x) = φ

tanh



x − x



where ξ

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

δφ

(

)), the

corresponding change in F is

δF =

δ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

δφ

(

) 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

had before, then we write

δF

δQ

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

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→ φ

= φ

(r) = φ(r − u).

Then

δφ(r) = φ

(r) − φ(r) = −u · ∇φ(r) + O(u

Then

δF =

δφ

δF

δφ

dr = −

µu · ∇φ dr

φ∇ · (µu) dr =

(φ∇µ) · u dr =

(φ∇

µ)u

dr.

using incompressibility.

We can think of the free energy change as the work done by stress,

δF =

(r)ε

(r) dr,

where

∇

is the strain tensor, and

is the stress tensor. So we can

write this as

δF =

∇

dr = −

(∇

dr.

So we can identify

∇

= −φ∇

µ.

So µ also contains the “mechanical information”.

3.2 Functional integrals

Given a coarse-grained ψ, we have can define the total partition function

−βF

TOT

= Z

TOT

−βF [ψ]

D[ψ],

where D[

] is the “sum over all field configurations”. In mean field theory, we

approximate this

TOT

by replacing the functional integral by the value of the

integrand at its maximum, i.e. taking the minimum value of

[

]. 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

]

of volume

with periodic boundary conditions. We can define the Fourier modes

√

ψ(r)e

−iq·r

dr,

Since we have periodic boundary conditions,

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”

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

= 0 for

|q| > q

max

. This leaves us with finitely many degrees of freedom.

The normalization of ψ

is chosen so that Parseval’s theorem holds:

|ψ|

dr =

|ψ

We can then define

D[ψ] =

dψ

Since we imposed a

max

, this is a finite product of measures, and is well-defined.

In some sense,

max

is arbitrary, but for most cases, it doesn’t really matter

what

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

, 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

’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

max

, and we take the product up to

there. In most of our calculations, we need such a

max

to make sense of our

integrals, and that will be left implicit. Most of the time, the results will be

independent of

max

(for example, it may give rise to a constant shift to

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

will still be complex in general, but they will

not be independent. Instead, we always have

= ψ

∗

−q

Thus, we should only multiply over half of the possible

’s, and we usually

denote this by something like

In practice, there is only one path integral we are able to compute, namely

when βF is a quadratic form, i.e.

βF =

φ(r)G(r − r

)φ(r

) dr dr

−

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 [ψ

] =

G(q)φ

−q

−

Example. We take Landau–Ginzburg theory and consider terms of the form

βF [φ] =

ξφ

− hφ +

(∇φ)

(∇

φ)

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 {φ

} =

(a + κq

+ γq

)φ

−q

−

(a + κq

+ γq

)φ

−q

−

So our G(q) is given by

G(q) = a + kq

+ γq

To actually perform the functional integral, first note that if

h 6

= 0, then we

can complete the square so that the

term goes away. So we may assume

= 0.

We then have

TOT

dφ

−βF {φ

}

dφ

−|φ

G(q)

Each individual integral can be evaluated as

dφ

−|φ

G(q)

ρ dρ dθ e

−G(q)ρ

G(q)

where φ

= ρe

iθ

. So we find that

TOT

G(q)

and so

βF

= −log Z

log

G(q)

We now take large V limit, and replace the sum of the integral. Then we get

βF

(2π)

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

i =

−k

−

−q

G(q)

dφ

We see that this is equal to

∂Z

∂G(k)

= −

∂ log Z

∂G(k)

G(k)

We could also have done this explicitly using the product expansion.

This

(

) is measured in scattering experiments. In our previous example,

for small k and κ > 0, we have

S(q) =

a + κk

+ γk

≈

−1

1 + k

, ξ =

where ξ is the correlation length. We can return to real space by

hφ

(r)i =



|φ(r)|



hφ

(2π)

max

a + κq

+ γq

4 The variational method

4.1 The variational method

The variational method is a method to estimate the partition function

−βF

TOT

−βF [φ]

D[φ]

when

is not Gaussian. To simplify notation, we will set

= 1. It is common

to make a notational change, where we replace

TOT

with

and

[

] with

[

We then want to estimate

−F

TOT

−F [φ]

D[φ].

We now make a notation change, where we write

TOT

, and

[

] as

[

]

instead, called the effective Hamiltonian. In this notation, we write

−F

−H[φ]

D[φ].

The idea of the variational method is to find some upper bounds on

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

[φ], and similarly define

−F

−H

[φ]

D[φ].

We can then write

−F

−H

D[φ]

−H

−(H−H

)

D[φ] = e

−F

−(H−H

)

where the subscript 0 denotes the average over the trial distribution. Taking the

logarithm, we end up with

F = F

− loghe

−(H−H

)

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

≥ hlog Y i

Applying this to our situation gives us an inequality

F ≤ F

− hH

− Hi

= F

− hH

+ hHi

= S

+ hHi

This is the Feynman–Bogoliubov inequality.

To use this, we have to choose the trial distribution

simple enough to

actually do calculations (i.e. Gaussian), but include variational parameters in

. We then minimize the quantity

− hH

hHi

over our variational

parameters, and this gives us an upper bound on

. We then take this to be

our best estimate of

. If we are brave, we can take this minimizing

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

. If we simply want the different parts to stay

away from each other, then we can have a configuration that looks like

In general, we expect that there is such an order along the

direction, as

indicated, while there is no restriction on the alignments in the other directions.

So the system is a lattice in the

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

for some

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

z for some A.

We again use Landau–Ginzburg model, which, in our old notation, has

βF =



(∇φ)

(∇

φ)



dr.

If we write this in Fourier space, then we get

βF =

(a + κq

+ γq

)φ

−q

−(q

)

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

= 0. Then this is a

Gaussian model with

G(q) = a + κq

+ γq

= τ + α(q − q

)

Varying

gives a linear shift in

(

). As we change

, we get multiple different

curves.

G(q)

a = a

Thus, as a decreases, S(q) = h|φ

i =

G(q)

blows up at some finite

q = q

−κ

2γ

, a

4γ

We should take this blowing up as saying that the

|q|

states are highly

desired, and this results in an ordered phase. Note that any

with

|q|

highly desired. When the system actually settles to an ordered state, it has to

pick one such

and let the system align in that direction. This is spontaneous

symmetry breaking.

It is convenient to complete to square, and expand

about

and

a = a

. Then we have

G(q) = τ + α(q − q

)

where

τ = a − a

, α =

) = −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

. As suggested before, we try a solution of the

form

φ = A cos q

which is smectic along

. We can then evaluate the free energy (per unit volume)

to be

βF

= βF [φ] =

(∇φ)

(∇

φ)

where the bar means we average over one period of the periodic structure. It is

an exercise to directly compute

, (∇φ)

, (∇

φ)

, φ

This gives

βF



+ κA

+ γA







(a − a

)

| {z }





Note that

is fixed by the system as we originally had, while

is the amplitude

of the fluctuation which we get to choose. Plotting this, we get a familiar graph

βF/V

τ > 0

τ < 0

τ > 0

τ < 0

τ >

0, then the optimum solution is given by

= 0. Otherwise, we should

pick

A 6

= 0. Observe that as we slowly reduce

across 0, the minimum varies

continuously with τ :

|A|

Variational method

We now consider the variational theory. In the notation we were using previously,

we have

H =

−q

G(q) +

−(q

)

Our trial H

−q

J(q).

Since this is a Gaussian model, we know that

log

J(q)

To use our inequality, we need to evaluate our other two bits. We have

hφ

−q

J(q).

We already calculated

hφ

−q

J(q)

Thus, we have

Here it is clear that we must impose a cutoff on

. We can think of this 1 as the

equipartition theorem.

We can also compute

hHi

J(q)

G(q) +

hφ

−(q

)

| {z }

In the Gaussian model, each

is a zero mean Gaussian random variables, which

have certain nice properties. Wick’s theorem tells us we have

habcdi

= habi

hcdi

+ haci

hbdi

+ hadi

hbci

Applying this, together with the result

hφ

= h|φ

,−q

we obtain

U = 3

h|φ

= 12

J(q)

Thus, we have

G(q)

J(q)

We can then compute

F =



log

J(q)

− 1 +

G(q)

J(q)



J(q)

We minimize over J(q) by solving

∂

∂J(q)

= 0

for all q. Differentiating, we obtain

J(q)

−

G(q)

J(q)

−

V J(q)

J(q

)

= 0.

Multiplying through by J

, we get

J(q) = G(q) +

J(q

)

For large V , we can replace the sum by an integral, and we have

J(q) = G(q) +

(2π)

J(q

)

It is very important that once we have fixed

(

), the second term is a constant.

Writing

C =

J(q

)

(2π)

J(q

)

we can then say the minimum is given by

(

) =

(

) + 3

. Thus, solving for

J(q) is equivalent to finding a C such that

C =

(2π)

G(q) + 3bC

This is a self-consistency equation for C (and hence J).

There is a very concrete interpretation of

. Recall that

hφ

is the average

value of φ

at some point r (which is independent of r. So we can compute

hφ(r)

i =

hφ(r)

i dr =

h|φ

i =

J(q)

Thus, what our computation above amounts to is that we replaced

(

) by

G(q) + 3bhφ

A good way to think about this is that we are approximating the

term in

the free energy by

≈

hφ

iφ

We have to pick a value so that this approximation is consistent. We view the

hφ

above as just another constant, so that the free energy is now a Gaussian.

We can compute the expectation of

using this free energy, and then we find

that

hφ

i =

(2π)

G(q) + 3bhφ

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

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

)

then we have

J(q) = ¯τ + α(q − q

)

, ¯τ = 3bhφ

The self-consistency equation now states

¯τ = τ +

(2π)

¯τ + α(q − q

)

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

= 3, we can

make the approximation.

(2π)

¯τ + α(q − q

)

2π

∞

¯τ + α(q − q

)

≈

3bq

2π

∞

¯τ + α(q − q

)

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

¯τ = τ +

√

¯τ

, s =

2π

√

∞

1 + y

∼

√

The precise value of

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|φ

i = S(q) =

¯τ + α(q − q

)

becomes large but finite. Note that since

¯τ

is finite, this sets some sort of length

scale where all q with |q − q

| ∼ ¯τ 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, φ = φ

+ δφ, where

= A cos q

We can do computations similar to what we did above by varying

δφ

, and obtain

(

). Then the global minimum over all

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)

For finite A, the quartic term now has an extra contribution, and we bet

¯τ = τ +

√

¯τ

3bA

Compare this with mean field theory, where we have

MF T

(A) =



τA



We see that for small

, the fluctuations are large, and mean field theory is

quite off. For large

, the fluctuation terms are irrelevant and mean field theory

is a good approximation. Thus, we get

MFT

We see that as

decreases, the minimum discontinuously jumps from

= 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|

We have not calculated

, but it can be done. We shall not do this, but

we can state the result (Brazovskii, 1975):

' −(sb)

2/3

, A

' s

1/3

−1/6

It turns out the variational approach finally breaks down for

τ < τ

∼ −s

3/4

1/2

We have

 τ

b 

√

. 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

is, we need to find leading corrections to

, as

Brazovskii did. In general, there is no reason to believe

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

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

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

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 =

L dt, L = T − V.

The equations of motion are

δA

δx(t)

∂L

∂x

−



∂L

∂x



= 0.

For example, if

L =

m ˙x

− V (x),

then the equation of motion is

−

δA

δx(t)

= m¨x + V

(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

= −ζ ˙x.

(ii) We introduce a noise f, with hfi = 0.

We set F

BATH

= F

+ 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

exp



−

f(t)



where

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

)i = σ

δ(t − t

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 (

, t

) and (

, t

). Given

any path

(

) 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

[x(t)] “=” P[f] = N

exp

−

2σ



ζ ˙x −

δA

δx



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

, denoting the forward probability. We can also consider

the backward probability

[

(

)], which is the probability of going along the

path in the opposite direction, from (x

, t

) to (x

, t

To calculate this, we use the assumption that

δA

δx

is time-reversal invariant,

whereas ˙x changes sign. So the backwards probability is

[x(t)] = N

exp

−

2σ



−ζ ˙x −

δA

δx



dt.

The point is that at equilibrium, the probability of seeing a particle travelling

along

(

) 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

[

(

)] is the same as

[

(

)]. Instead, if at

equilibrium, there are a lot of particles at

, then it should be much less likely

for a particle to go from x

to x

than the other way round.

Thus, what we require is that

[x(t)]e

−βH

= P

[x(t)]e

−βH

where

(

, ˙x

). 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

[x(t)]

exp



−

2σ



(ζ ˙x)

− 2ζ ˙x

δA

δx



δA

δx





exp



−

2σ



(ζ ˙x)

+ 2ζ ˙x

δA

δx



δA

δx





= exp



−

2σ



−4ζ ˙x

δA

δx



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 =

˙x

+ V (x).

We then find that



˙x

∂L

∂ ˙x



−

= ¨x

∂L

∂ ˙x

+ ˙x

∂L

∂ ˙x

−



˙x

∂L

∂x

+ ¨x

∂L

∂ ˙x



= ˙x



∂L

∂ ˙x



−

∂L

∂x



= −˙x

δA

δx

Therefore we get

˙x

δA

δx(t)

dt = −(H(x

, ˙x

) − H(x

, ˙x

)) = −(H

− H

Therefore, the principle of detailed balance tells us we should pick

= 2k

T ζ.

This is the simplest instance of the fluctuation dissipation theorem.

Given this, we usually write

f =

T ζΛ,

where Λ is a Gaussian process and

hΛ(t)Λ(t

)i = δ(t − t

We call Λ a unit white noise.

In summary, for a particle under a potential V , we have an equation

m¨x + V

(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

)i = σ

δ(t − t

) = 2k

T ζδ(t − t

In general, in the coarse grained world, suppose we have mesostates

A, B

with probabilities e

−βF

and e

−βF

, then we have an identical statement

−βF

P(A → B) = e

−βF

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

, how does

this evolve over time?

It is convenient to consider the overdamped limit, where

= 0. Our equation

then becomes

ζ ˙x = −∇V +

T ζΛ.

Dividing by ζ and setting

M = ζ

−1

, we get

˙x = −

M∇V +

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 ∆

in a time interval

∆t. Equivalently, we are asking Λ to change by

∆Λ =

√

T ζ

(ζ∆x + ∇V · ∆t).

Thus, the probability of this happening is

W (∆x, x) ≡ P

∆t

(∆x) = N exp



−

4ζk

T ∆t

(ζ∆x + ∇V ∆t)



We will write

for ∆

. Note that

(

u, x

) is just a normal, finite-dimensional

Gaussian distribution in

. We can then calculate that after time ∆

, the

expectation and variance of u are

hui = −

∇V

∆t, hu

i − hui

∆t + O(∆t

We can find a deterministic equation for P (x, t), given in integral form by

P (x, t + ∆t) =

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 − u∇P +

∇



W (u, x) − u∇W +

∇



where all the gradients act on

, not

. Applying the integration to the expanded

equation, we get

P (x, t + ∆t) = P (x, t) − hui∇P +

i∇

P − P ∇hui,

Substituting in our computations for hui and hu

i gives

P (x, t)∆t =



∇V

∇P +

∇

P +

P ∇



∆t.

Dividing by ∆t, we get

P =

∇

P +

∇(P ∇V )

= D∇

P +

M∇(P ∇V ),

where

D =

M =

are the diffusivity and the mobility respectively.

Putting a subscript

to emphasize that we are working with one particle,

the structure of this is

= −∇ · J

= −P

D∇(log P

+ βV )

= −P

M∇µ(x),

where

µ = k

T log P

+ V

is the chemical potential of a particle in V (x), as promised. Observe that

– This is deterministic for P

–

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

(

) as

long as the temperature is constant, i.e.

(

) =

D(x)

. In this case,

the Langevin equation says

˙x = −

M(x)∇V +

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

non-interacting colloids under the same potential

(

). 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

, h˙ρi = N

Then our previous discussion entails ˙ρ evolves by

˙ρ = −∇ · J,

where

hJi = −ρ

M∇µ, µ = k

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

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

is a noise current. This is the Langevin equation for a fluctuating field

ρ(r, t).

We can fix the distribution of

by requiring detailed balance as before. We

will implement this for a constant

, 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



−

2σ

dr |j(r, t)|



This corresponds to

(r, t)j

, t

)i = σ

δ(r − r

)δ(t − t

We now repeat the detailed balance argument to find σ

. We start from

J + M∇µ = j.

Using F to mean forward path, we have

[J(r, t)] = N exp



−

2σ

dr |J + M∇µ|



where

µ =

δF [ρ]

δρ

We consider the backwards part and get

[J(r, t)] = N exp



−

2σ

dr | − J + M ∇µ|



Then

log

= −

dr J · ∇µ.

We integrate by parts in space to write

dr J · ∇µ = −

dr (∇ · J)µ =



˙ρ

δF

δρ



dF [ρ]

So we get

log

= −

− F

So we need

= β,

or equivalently

= 2k

T M.

So our final many-body Langevin equation is

˙ρ = −∇ · J

J = −M ∇



δF

δρ



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 = −

M∇V +

MΛ,

and we turned this into a Fokker–Planck equation

P = −∇· J

J = −P

M∇µ

µ = k

T log P + V (x).

We then write this as

P = ∇·



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

(

) with

(

r, t

), and we replace

(

x, t

)

with

[

(

);

]. We then replace

∇

with

δρ(r)

. So the Fokker–Planck equation

becomes

P [ρ(t); t] =

δρ



T ∇ ·

M∇



δρ

+ β

δF

δρ



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

(

), and

evolves under the equations

φ = −∇ · J

J = −M ∇µ +

T MΛ

µ =

δF

δφ

Recall that the system is modelled under the Landau–Ginzburg free energy

F [φ] =



| {z }

(∇φ)



dr.

We then have

µ = aφ + bφ

− κ∇

φ.

As before, the mean field theory for equilibrium looks like

a > 0

a < 0

a > 0

a < 0

Here

±φ

are the spinodals, where the second derivative changes sign. This

gives the phase diagram

−1 1

a(T ) = 0

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

φ| < φ

, and

(

)

0 implies local instability. The

system rapidly becomes phase separated. This is spinodal behaviour.

–

In (3), we have

< |

φ| < φ

. 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

∂φ

− κ∇

φ = f

(

φ) +

φf

(

φ) − κ∇

φ.

Note that the first term is a constant. We then have

φ = −∇ · J

J = −M ∇[f

φ − κ∇

φ] +

T MΛ.

We drop the tildes and take the Fourier transform to get

= −Mq

+ κq

)φ

+ iq ·

T MΛ

Compare this with an overdamped particle in a simple harmonic oscillator,

V =

κx

where we have

˙x = −

Mκx +

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φ

(0)φ

−q

(t)i = S(q)e

−r(q)t

This

(

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

= −Mq

+ κq

)

| {z }

r(q)

+ iq ·

T MΛ

but crucially, now

(

)

0, so it is possible to have

(

)

0. The system is

unstable.

If we ignore the noise by averaging the equation, then we get

i = −r(q)hφ

So if we have a noisy initial state φ

(0), then the perturbation grows as

hφ

(t)i = φ

(0)e

−r(q)t

When

(

)

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:

r(q)

∗

The maximally unstable mode

∗

is given by the solution to

(

∗

) = 0, which

we can easily check to be given by

∗

−f

2κ

Now consider the equal time non-equilibrium structure factor

(t) = hφ

(t)φ

−q

(t)i ∼ S

(0)e

−2r(q)t

As time evolves, this gets more and more peaked around q = q

∗

(t)

∗

So we see a growth of random structure with scale

L ∼ π/q

∗

. This process is

called spinodal decomposition.

Note that this computation was done on the assumption that

is small, where

we ignored the quadratic terms. At intermediate

, 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

, where

= f

(2π)

max

(t) d

This says

is less negative as the fluctuations grow. Since

∗

−

this moves to a smaller

. So

(

)

∼ π/q

∗

(

) starts increasing. This is called

domain growth.

In the late stages, we have large regions of

φ ≈ ±φ

, 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

σA(t)

where

(

) 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



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) ∼





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

−φ

, where

small. The system is locally stable, so

(

)

0 for all

. 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



−

|J + M∇µ|

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

−φ

then at some point, there will be some random small droplet with

= +

small radius

. This is already unlikely, but after this, we need

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

. First there is

the cost of having an interface, which is 4

πσR

. However, the fact that we have

areas is energetically favorable, and this grows as the volume. So we get a

cubic term ∼ −R

. If we add these two together, we get a barrier:

F (R)

∗

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

±φ

, 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

−φ

Then in the phase separated state, the +

phase will lump together in a droplet

(if we had

φ = 0, then we would have a horizontal interface).

Within each region 1 and 2, the value of

is constant, so the term that contributes

to the free energy is

f(φ) =

We would expect 1 and 2 to respectively be located at

1 2

When we have a spherical interface, 1 and 2 are not exactly at

±φ

. 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

= µ

bulk

, Π

bulk

= Π

bulk

This means the points 1, 2 have a common tangent

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 (Π

−

)

πR

while the interfacial tension pulls the boundary down by 2

πRσ

. In general, in

dimensions, we require

= Π

(d − 1)

This is called the Laplace pressure.

In static equilibrium, we still require

, since this is the same as saying

∇µ

= 0. So

has the same slope at

and

. However, the two tangent

lines no longer have a common intercept, but they are separated by

(

d −

1).

So it looks like

−Π

To solve for this, we take the approximation that

is small for

decently

large. Then we can write

= f(−φ

+ δ

) ≈

(−φ

)δ

+ f (−φ

)

= f(+φ

+ δ

) ≈

(+φ

)δ

+ f (+φ

So µ

= αδ

, where α = f

(±φ

). So we find that up to first order, δ

= δ

To compute δ, we compute

= µ

− f

= −αδφ

Similarly, we have Π

= +αδφ

. So

− Π

= −2αφ

δ.

Since this equals −(d − 1)

, we have

δ =

d − 1

2αφ

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

inside a bath with

−φ

where

ε 6

(

). This

is called the supersaturation. Note that to have a

droplet of radius

, the value of

inside and immediately outside the droplet

must be

±φ

. Outside of the droplet, the value of

will slowly decay to

−φ

+ ε. Thus, outside of the droplet, we write

φ(r) = −φ

φ(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

(

) at any point in time, we set

= 0.

So ∇

φ = 0. We solve this with boundary conditions

φ(∞) = ε,

φ(R

) = δ.

So we have

φ = ε + (δ − ε)

Now if we assume this is what

(

) looks like, then the current just outside the

droplet gives

J(R

) = −M∇µ = −αM

∂

∂r



= αM(δ − ε)



r=R

αM(δ − ε)

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

. So mass

conservation implies

2φ

R = −J = −

αM(δ − ε)

Thus, we conclude that

R =

2φ



αM

(ε − δ(R))



We can plug in our previous expression for

. Fixing

, we can plot

as follows:

∗

where

∗

αεφ

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 ≈

2φ

αM

(ε − δ(

R)).

We know that the value of

is also determined by

, so we know

ε − δ

(

) is of

order δ(

R). Hence

R ∼

Mσ

∼

Mσt

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σ



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

particles inside the droplet exerting ideal gas pressure, then we

have

− Π

2σ

4/3πR

We again have µ

= µ

. This ends up giving a new boundary condition at R

φ(R

) =

αRφ

−

3Nk

8αφ

πR

2αφ

(Π

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 =

2φ



αM



ε −

αφ

3Nk

8αφ

πR



We now see that there is a stable fixed point

, 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

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

, then our

equation is now

φ + v · ∇φ = −∇ · J.

The

v · ∇φ

is called the advection term. Our current

is the same as before,

with

J = −M

δF

δφ

T MΛ.

We also need an evolution equation for

, 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

= Σ

+ Σ

– The Σ

term is the pressure term, given by

= −P δ

We should think of this P as a Lagrange multiplier for incompressibility.

– The Σ

term is the viscous stress, which we can write as

= η(∇

+ ∇

)

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

= −Πδ

− κ(∇

φ)(∇

φ), Π = φµ − 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Σ

(r, t)Σ

, t

)i = 2k

T η



+ δ

−



δ(r−r

)δ(t−t

The last term

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 + η∇

v + −φ∇µ + ∇ · Σ

Hence Model H has equations

φ + v · ∇φ = −∇ · J

J = −M ∇µ +

T MΛ

∇ · v = 0

ρ(

v + v · ∇v) = η∇

v − ∇P − φ∇µ + ∇ · Σ

Compared to Model B, we have the following new features:

(i) −φ∇µ drives deterministic fluid flow.

(ii) Σ

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

= 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

i ∼ Dt

, where

D =

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

(

), which determines the size and separation of

droplets. We can then calculate the collision time

∆t '

∼

Each collision doubles the volume, and so

R →

1/3

. Taking the logarithm, we

have

∆ log

R ∼

log 2

in time ∆t.

So we crudely have

∆ log

∆t

∼

log 2

If we read this as a differential equation, then we get

d log

R ∼

So we find that

R ∼

R(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σ



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.

Again assume we have a single length scale

(

), given by the domain size.

Again we assume we have a single length scale. As time goes on, we expect

(

)

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

. Since there is only one

length scale, we must have

K ∼ v.

The Laplace pressure then scales as ∼

The noise terms Σ

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) = η∇

v − ∇P − φ∇µ,

We make single length scale approximations as before, so that

v ∼

and

∇ ∼ L

−1

. Then we have an equation of the form.

L + ρ

∼ η

+ Lagrange multiplier +

(∗)

where we recall that at curved interfaces,

µ ∼ ±

. Here we have a single variable

(

), and three dimensionful parameters

ρ, η, σ

. We can do some dimensional

analysis. In d dimensions, we have

ρ = M L

−d

, η = ML

2−d

−1

, σ = ML

3−d

−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

ρσ

, t

ρσ

One can check that these are the only combinations with units

L, T

. So we must

have

L(t)

= f





We now substitute this into the equation (

∗

). We then get a non-dimensionalized

equation for (∗)

αf

+ βf

/f = γ

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

(or small f ). Then we get

γf

= 0,

is a constant, and so

grows linearly with

. Putting all the

appropriate constants in, we get

L ∝

This is called the viscous hydrodynamic regime, VH .

(ii)

For large

, we assume we have a power law

(

)

∼ x

, where

y >

0 (or

else f would not be large). Then

¯αx

y−2

βx

y−2

= ¯γx

−y−1

+ δx

−2y

It turns out at large

, the

−y−1

term is negligible, scaling wise. So we

have y − 2 = −2y, or equivalently, y =

. So

∼





2/3

Putting back our factors, we have

L ∼





1/3

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

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

ψ ∼

6, then

we always get a bicontinuous medium. If

ψ <

3 or

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

is strictly symmetric (

(

) =

(

−φ

)), the

only case where we have bicontinuous phase is on the ψ =

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 =



|p|

(∇

)(∇

)



dr ≡

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

|∇ · p|

p · ∇ ∧ p|

p ∧ p|

Here

is not conserved, so

is not of the form

−∇ · J

. Instead, (without flow)

we have

p = −Γh, h =

δF

δp(r)

where

(

) 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

= −Γh,

where D is some sort of comoving derivative. For the scalar field, we had

Dφ

φ + v · ∇φ.

Note that the advective term is trilinear, being first order in

∇

and

. For

there is certainly something like this going on. If we have a translation, then

gets incremented by

∆p = v · ∇p ∆t,

as for a scalar.

There is at least one more thing we should think about, namely if

rotational, then we would expect this to rotate

as well. We have the corotational

term

∆p = ω ∧ p ∆t,

where ω is the angular velocity of the fluid, given by

ijk

Ω

, Ω

(∇

− ∇

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

∆p = −ξD · p∆t

with D

(

∇

) and

a parameter. This is the irrotational part of

the derivative. The reason is that in a general flow,

needn’t simply rotate with

ω. Instead, it typically aligns along streamlines. In total, we have

= (∂

+ 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

= −Γh, h =

δF

δp

We next need an equation of motion for v. We can simply

ρ(∂

+ v · ∇)v = η∇

v − ∇P + ∇ · Σ

where Σ

is a stress tensor coming from the order parameter.

To figure out what Σ

should be, we can consider what happens when we

have an “advective” elastic distortion. In this case, we have

= 0, so we have

p = −v · ∇p − Ω · p + ξD · p,

The free energy change is then

δF =

δF

δp

p ∆t dr = ∆t

h · p dr,

On the other hand, the free energy change must also be given by

δF =

∇

(r) dr,

the product of the stress and strain tensors. By matching these two expressions,

we see that we can write

= Σ

(1)

+ Σ

(2)

+ Σ

(3)

where

∇

(1)

= −p

∇

, Σ

(2)

− p

), Σ

(2)

+ p

Analogous results for a nematic liquid crystal holds, whose derivation is a

significant pain. We have

= −ΓH, H

δF

δQ

−



δF

δQ



The second term in H is required to ensure Q remains traceless.

The total derivative is given by

= (∂

+v ·∇)Q+ (Ω ·Q −Q ·ω) + ξ(D ·Q + Q ·D)−2ξ



Q +



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,1

+ Σ

Q,2

where

∇

Q,1

k,`

= −Q

∇

Q,2

= Q · H − H · Q − ξ



H + 2

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

= 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

(

). If the quench is deep enough, we have spinodal-like

instability, and we quickly get locally ordered. Coordinate-independently, we

can write

Q =



−



Since all this ordering is done locally, the principal axis

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

±φ

, but we get

domain walls between the ±φ

phases.

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 = −

q = +

The charge

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

−

, the order

parameter made a half turn in the opposite sense to how we moved around the

defect. In the q = +

case, they turned in the same sense.

We see that

q ±

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

= +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

, then they tend to dissociate

into multiple

defects. This is due to energetic reasons. The elastic

energy is given by

ell

|∇ · Q|

∼

|(∇ · n)n + n · ∇n|

If we double the charge, we double the

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).

This

(

) is random, and tend to vary continuously. However, topological defects

are present, which cannot be ironed out locally. All topological defects with

|q| >

dissociate quickly, and we are left with q = ±

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

˜κ

|(∇ · n)n + ∇n|

where ˜κ = κλ

. Dimensionally, we have

∇ ∼

So we have an energy

E ∼ ˜κ

dr ' ˜κ log





where

is the mean spacing and

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

∝

˜κ

. Under this force, the

defects move with overdamped motion, with the velocity being proportional to

the force. So

R ∼

L ∝

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

˜κ

and so and so we similarly have

L(t) ∼ t

1/2

However, in three dimensions, the

defects are the same topologically. In

other words, we can change +

−q

via continuous, local deformations. This

involves rotating out to the

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

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

dimensions, we have

d−1

, the (d − 1)-dimensional unit sphere.

Example.

For nematic liquid crystals in

-dimensions, we have

d−1

which is obtained from S

d−1

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

that encloses a defect core, and consider the map

D → M

that assigns to each point the order parameter at that point. In our cases,

a circle S

, and so f is a loop in M.

We say two mappings

, f

, are homotopic if they can be continuously

deformed into each other. Defects lie in the same homotopy class if maps for

all

’s enclosing them are homotopic. The fundamental group

(

) is the set

of all homotopy classes of maps

→ 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

→ 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 π

n−1

(M).

Example.

Observe that

is actually just

in disguise, and so

(

) =

The generator of π

(RP

) is the charge

topological defect.

Example. We can visualize RP

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

, and this is the q =

defect.

Observe that in the two-dimensional case, the

defects correspond to

going along the top arc and bottom arc from the left to right respectively. In

, there is then a homotopy between these two paths by going through the

disk. So in RP

, they lie in the same homotopy class.

In general, it is easy to see that

(

) =

, so

is the unique

non-trivial defect.

This is particularly interesting, because two

defects can merge and

disappear! Similarly, what you would expect to be a

= 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

(

)

∼

. 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

. Platinum is a catalyst of the decompo-

sition

→ 2H

O + O

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 =



(∇φ)



dr.

To model our system without time reversal symmetry, we put

˜µ =

δF

δφ

+ λ(∇φ)

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)

(

∇φ

)

is not the functional derivative of any

. This breaks the free energy

structure, and

6= e

−β(F

−F

)

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(φ) =

δφ

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)

= (µφ − f)

+ ∆.

This is found by solving J = 0, so

˜µ =

∂f

∂φ

− κ∇

φ + λ(∇φ)

= const.

(vi) There is a further extension, active model B+, where we put

J = −∇˜µ +

√

2DΛ + ζ(∇

φ)∇φ.

This extra term is similar to

∇

(

∇φ

)

) 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

. 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 =



(∇

)(∇

)



dr.

As for active model B,

can acquire gradient terms that are incompatible with

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

, then it is natural that

wants to

translate along its own direction at some speed

. Thus,

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

= −Γh,

where



∂

∂t

+ v · ∇



p + Ω · p − ξD · p + wp · ∇p.

The Navier–Stokes/Cauchy equation is now

(∂

+ v · ∇)v = η∇

v − ∇P + ∇ · Σ

(p)

+ ∇ · Σ

where as before,

∇ · Σ

(p)

= −p

∇

+ ∇



− p

) +

+ p

)



and we have a new term Σ

given by the active stress, and the lowest order

term is

ζp

. This is a new mechanical term that is incompatible with

. We

then have

∇ · Σ

= (∇ · 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

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

= −ΓH, H =



δF

δQ



traceless

where

is given by

= (∂

+ v · ∇)Q + S(Q, K, ξ).

Here K = ∇v and

S = (−Ω · Q − Q · Ω) − ξ(D · Q + Q · D) + 2ξ



Q +



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

= ζ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 = ∇ · Σ

= ζ∇ · 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

= yg. Then the viscous stress is

= ηg



0 1

1 0



We have

∝ ζλ



nn −



= ζλ



0 1

1 0



is at 45

◦

exactly. Note that the sign of

affects whether it reinforces or

weakens the stress.

A crucial property is that Σ

does not depend on the shear rate. So in the

contractile case, the total stress looks like

TOT

In the extensile case, however, we have

TOT

∗

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 = −

q = +

Note that the

−

charge is symmetric, and so by symmetry, there cannot be

a net body stress. However, in the

= +

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

are formed locally. The

−

stay put

and the +

ones self-propel, and depending on exactly how the defect pairs are

formed, the +

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).