VOLUME 87, NUMBER 7

PHYSICAL REVIEW LETTERS

13 AUGUST 2001

Simulating Particle Dispersions in Nematic Liquid-Crystal Solvents
Ryoichi Yamamoto*Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, United Kingdom(Received 15 April 2001; published 27 July 2001)
A new method is presented for mesoscopic simulations of particle dispersions in nematic liquid-crystalsolvents. It allows efficient first-principles simulations of the dispersions involving many particles withmany-body interactions mediated by the solvents. A simple demonstration is shown for the aggregation
process of a two dimensional dispersion.
DOI: 10.1103/PhysRevLett.87.075502

PACS numbers: 61.30.Cz, 61.20.Ja, 61.30.Jf

Dispersions of small particles in host fluids such as colloidal suspensions and emulsions are of considerable technological importance, and often appear in our everyday
life in paints, foods, and drugs. Many kinds of exotic interactions between particles mediated by the host fluids arepossible, including screened Coulombic [1], depletion [1],fluctuation induced [2], and surface induced [3] forces. Astriking example occurs when spherical particles are immersed in a liquid-crystal solvent in the nematic phase. Fora single particle, the orientation of the solvent molecules is
distorted due to the anchoring of the solvent molecules atthe particle surface. Extensive studies have been done on
this effect, and several characteristic configurations of thenematic field around a spherical particle have been identified [4–9]. When the strength of anchoring is increasedso that normal anchoring is preferred, the solvent changes
from quadrupolar to dipolar symmetries around the particle. When more than two particles are immersed in thesolvent, long-range anisotropic interactions are induced between particles due to elastic deformations of the nematicfield [10–13]. The anisotropic interactions can have a pronounced effect not only on the local correlations of the particles [10], but also on their phase behavior [14–17] andon their mechanical properties [15].
Since analytical approaches for investigating these kindsof complex materials are extremely difficult, computer
simulations are essential to investigate their static and dynamical properties. In most dispersions, the host fluid
molecules are much smaller and move much faster thanthe dispersed particles. This enables us to assume that thehost fluid is in local equilibrium for any given particle configurations, and thus it is usually a good idea to use somecoarse grained mesoscopic descriptions for the host fluids
rather than treating them as fully microscopic molecules[18]. In the case of charged colloidal suspensions, a mesoscopic method for the first-principles simulations can bederived by treating the counterions as a charge density[19]. For the particle dispersions in nematic solvents considered here, the mesoscopic coarse grained free energy
for the nematic solvent is well known to be the Frank freeenergy [20], and the total free energy F of the system consists of the following two parts: the bulk term Fel which
presents elastic energy of the nematic and the surface term

Fs which determines anchoring of the nematic field at theparticle surface. Let n共r兲 be the director, a common direction on which solvent molecules are aligned on averagewith a constraint jn共r兲j 苷 1. F can be given by functionals of 关n共r兲兴 for a given particle configuration 兵R1 · · · RN 其:KZF 共关n共r兲兴; 兵R1 · · · RN 其兲 苷
dr 关共= ? n兲2 1 共= 3 n兲2 兴
2WIdS 关1 2 共n ? n兲2 兴 , (1)12075502-1© 2001 The American Physical Society

0031-9007兾01兾 87(7)兾075502(4)$15.00

where K is the Frank constant with the single elastic constant approximation, W is the surface anchoring constant,and n is the unit vector normal to the colloidal surface[6,7]. The saddle-splay elastic term [7] is not consideredin Eq. (1). The integral in the first term, Fel , runs over
the whole solvent volume excluding the particles, and thatin the second term, Fs , runs over all solvent-particle interfaces. A simple scaling argument tells us Fel ~ Kad22and Fs ~ Wad21 with a and d being the particle radius
and the system dimension, respectively, thus the physics
should be determined by the ratio Fs 兾Fel ~ Wa兾K.
Although this type of free energy functional is sufficientfor performing Monte Carlo simulations where onlyvalues of F are needed for a given particle configuration, it is not useful for molecular dynamics (MD) or
Brownian-type simulations because the coupling between
solvent and the particles is given implicitly by limiting
the integration space in both Fel and Fs . This produces
mathematical singularities at the interface when one
calculates the force, fiPS 苷 2≠F 兾≠Ri , acting on eachparticle mediated by the nematic solvents. Calculating
the force is crucial for performing efficient simulationsof many particle systems. Another serious problem ofthis type of functional is that in order to give correctboundary conditions at the particle-solvent interface,
one has to use appropriate coordinates for performing
grid-based numerical simulations rather than the usual
Cartesian coordinates. This is generally difficult for particles with nonspherical shapes or for systems involving
many particles even when each particle has a sphericalshape. Also this makes the use of the periodic boundary
condition difficult.
075502-1

VOLUME 87, NUMBER 7

PHYSICAL REVIEW LETTERS

13 AUGUST 2001

To overcome these problems, we have modified Eq. (1) by using a smooth interface between the solvent and the particles so that the coupling is given explicitly in the integrants through the interface. The new free energy functional wepropose is!
√NXK ZF 共关q共r兲兴; 兵R1 · · · RN 其兲 苷
dr 1 2fi 共r兲 tanh关Rc2 共=a qbg 兲2 兴
4Rc2i苷1
∏N ∑X
Wj Z
d21
共=a fi 兲2 2 共=a fi 兲 共=b fi 兲qab ,(2)1dr2di苷1
This form is very convenient because we can computewhere a, b, g [ x, y, z and the summation conventionboth ≠fi 兾≠Ri and ≠共=a fi 兲兾≠Ri at any time since fi
is used. The explicit form of the interfacial profile fi
between dispersed particles and solvents is given byis an analytical function of Ri . (iii) Finally, we updateµ∂
the particle positions according to appropriate equations
a 2 jr 2 Ri j
1of motion such astanh
11 ,fi 共r兲 苷(3)2j
d 2 Ri苷 fiPP 1 fiPS 1 fiH 1 fiR ,(7)miwith the particle radius a and the interface thickness
dt 2j. Note that this reduces to Eq. (1) if Rc , j ! 0. Very
where fiPP is the force due to direct particle-particle interacrecently, a similar idea of using smooth interface was protions (hard or soft sphere for instance), and fiH and fiR are
posed for treating the hydrodynamic forces acting on parthe hydrodynamic and random forces. Repeating the steps
ticles dispersed in simple liquids [21]. In our case, the free(i)–(iii) enables us to perform first-principles mesoscopic
energy is given by functionals of a traceless and symmetsimulations for the dispersions containing many particlesric second-rank tensor qab 共r兲 苷 na 共r兲nb 共r兲 2 dab 兾dwithout neglecting many-body interactions.
rather than the director n共r兲 to take into account the symWe have performed simple simulations for a two dimenmetry of the nematic director 1n $ 2n automatically.
sional(2D) system to demonstrate our simulation proce2
2
The semiempirical functional form 1兾Rc tanh关Rc · · ·兴 isdure.
Thesystem has 100 3 100 lattice sites in a squareapplied in Eq. (2) to avoid the mathematical divergenceboxwithalinear length L 苷 100. Other physical parameof the elastic free energy density at the defect centers andtersarechosen
rather arbitrarily as Rc 苷 1, a 苷 5, and2
to limit its value to Df ⬃ K兾Rc , which is the correct
j苷
2,
where
theunit of length is the lattice spacing l.
energy density difference between isotropic and nematic
Since
thenematicconfigurationsin 2D can be expressedstates, in the defect core regions of size Rc . Another waybya
singlescalarfield
关u共r兲兴,thetilt angle of the director
to avoid the divergence would be to use the Landau –against
thehorizontal
共x兲
direction,
Eq. (4) then reduces tode Gennes–type free energy with an order parameter
≠q共r兲
dFdF
abQab 共r兲 苷 Q共r兲qab 共r兲, but this requires a prohibitively
苷苷 0,(8)small lattice spacing near the defect points [22].
du共r兲
≠u共r兲 dqab 共r兲
The simulation procedure is as follows. (i) For a given
with qxx 苷 cos2 u 2 1兾2, qyy 苷 sin2 u 2 1兾2, and qxy 苷particle configuration 兵R1 · · · RN 其, we obtain the interface
qyx 苷 cosu sinu. The boundary condition is fixed atprofile fi 共r兲 by Eq. (3). Then we can calculate the stableu共r兲 苷 0 at the edge of the box to avoid rotations of the共0兲
(or metastable) nematic configurations 关qab 共r兲兴 which satreference frame. We first calculated stable nematic conisfy the equilibrium condition
figurations around a single particle for different Wa兾K,Çand found two stable configurations. The first configuradF
苷0(4)tion, which we refer to as weak anchoring, contains no共0兲dqab 共r兲 关qab 共r兲兴苷关qab 共r兲兴
topological defect. In the second configuration, which weunder the director constraint 关na 共r兲兴2 苷 1. One canrefer to as strong anchoring, the particle is accompanied
perform this by numerical iterations such as the steepest
by two 21兾2 charge point defects. Typical examples ofdescent or the conjugate gradient method. (ii) Once
the weak and strong anchoring are shown in Fig. 1(a) with共0兲
Wa兾K 苷 10 and Fig. 1(b) with Wa兾K 苷 20, respec关qab 共r兲兴 is obtained, the force acting on each particle
tively. The distance between the defects and the particle
mediated by the nematic solvents follows directly from
center is about 1.3a. We note that both configurations
the Hellmann-Feynman theorem,
possess quadrupolar symmetries, and the latter would
共0兲
correspond to the Saturn ring configuration in three共r兲兴;兵R···R其兲≠F共关q1
Nab(5)
fiPS 共兵R1 · · · RN 其兲 苷 2
dimensional (3D) systems. Although in principle particles≠Rican be accompanied by one 21 charge hedgehog defectin 2D as well as in 3D, such configurations are unstable in≠fiK Z共0兲dr
苷tanh关Rc2 共=a qbg 兲2 兴
the present 2D system since the elastic penalty of having
24Rc≠Rm point defects with charge c scales as mKc2 . This wasZ
≠共=a fi 兲共0兲
directly confirmed by recent simulations with perfect
共=b fi 兲qab . (6)1 W j dr≠Rinormal anchoring [22] and also by our simulations. The
075502-2

075502-2

VOLUME 87, NUMBER 7(a)when the solvent is in the isotropic phase 共K 苷 0兲. We
then set K 苷 1 and calculated fiPS according to the presentprocedure. The particle configurations were updated bynumerically solving the steepest descent equation,(b)zFIG. 1. Director configurations around a single particlefor (a) the weak anchoring case without defect obtained atWa兾K 苷 10 and (b) the strong anchoring case accompaniedby two 21兾2 charge point defects indicated with the crosses
obtained at Wa兾K 苷 20. The white disks indicate the particles
with radius a 苷 5. Only 9% of the total system is shown for
display purpose.

total free energies F 兾W are plotted in Fig. 2 as functionsof Wa兾K.
While both configurations can coexist in
the narrow transition regime 11 # Wa兾K # 12.5, ourmodel predicts a clear first-order transition from the weakanchoring to the strong anchoring around Wa兾K ⯝ 11.5.
We next simulated the aggregation and ordering processof 30 colloidal particles after the isotropic to nematic transition of the solvent occurred. Here we used the periodic
boundary condition and set Wa兾K 苷 20 so that each particle is accompanied by two 21兾2 charge defects. Otherparameters are the same as in the previous single particle
case. The simulation was performed starting from a random particle configuration which is a typical configuration

3transition region
2F/WWeakanchoring

Strong anchoring

1010

1102

Wa / K
FIG. 2. The total free energies for the single particle cases
as functions of the strength of the anchoring constant Wa兾K.
Our model predicts a first-order transition from weak to strong
anchoring around Wa兾K ⯝ 11.5.

075502-3

13 AUGUST 2001

PHYSICAL REVIEW LETTERSdRi
苷 fiPS 1 fiPP ,dt(9)which is obtained by simply substituting d 2 Ri 兾dt 2 苷 0,fiR 苷 0, and fiH 苷 2z dRi 兾dt in Eq. (7). z 苷 1 is a friction constant and thus the off-diagonal components of
the hydrodynamic interaction was not considered. Here
we obtain fiPP 苷 2≠EPP 兾≠Ri from the repulsive
PN21 PNpart of
the Lennard-Jones potential, EPP 苷 0.4 i苷1 j苷i11 3
关共2a兾jri 2 rj j兲12 2 共2a兾jri 2 rj j兲6 1 1兾4兴 truncated at
the minimum distance jri 2 rj j 苷 27兾6 a, to avoid the particles overlapping each other within the core radius ⯝a. Snapshots from the present simulation are shown inFig. 3(a) for an aggregation stage, and in Fig. 3(b) at alater time, where the particles are forming ordered clusters
due to the anisotropic attractions between them. Note isadded that only up to two particle simulations have been
done so far [7] and simulations of more than three particles
would be extremely difficult or almost impossible by othermethods ever proposed.
In summary, we have developed an extremely powerful
simulation method to investigate particle dispersions interacting via anisotropic solvents. We proposed a free energyfunctional which is suitable for MD-type simulations. The
following modifications have been made to the usual Frankfree energy functional. (i) The free energy is given by
a functional of a tensor q rather than a vector n to take
into account symmetry of the nematic director 1n $ 2n.
(ii) The coupling between the nematic solvent and particles at the interfaces is introduced explicitly through asmooth interface so that we can analytically calculate theforce acting on each particle mediated by the host by taking derivatives of the free energy according to the particle positions. (iii) The value of the free energy density
is limited semiempirically to avoid a mathematical divergence in the defect centers. We have performed demonstrations for a 2D dispersion and confirmed that the method
works well even when the system contains point defects.
Applications of this method to 3D systems should haveno theoretical difficulties, but require heavier computation. This should allow the simulation of the chaining ofthe particles caused by the possible dipolar symmetry ofthe nematic configurations around a single particle. Although we have shown only simple demonstrations of themethod by performing simulations of the 2D system in thisLetter, simulations with physically more interesting situations such as systems with noncircular particles, asymmetric particle pairs with different particle size, or particles
with non-normal anchoring as well as more realistic simulations in 3D systems are now underway.
075502-3

VOLUME 87, NUMBER 7

PHYSICAL REVIEW LETTERS

13 AUGUST 2001

The author thanks Professor J. P. Hansen and Dr. A.
Louis for helpful discussions. He acknowledges also theMinistry of Education, Culture, Sports, Science and Technology of Japan for supporting his stay in Cambridge in2000兾2001. Calculations have been carried out at the Human Genome Center, Institute of Medical Science, University of Tokyo, and the Supercomputer Center, Institute ofSolid State Physics, University of Tokyo.

FIG. 3. The aggregation and ordering process of colloidal particles when the solvent exhibits the isotropic 共K 苷 0兲 to nematic
共K 苷 1兲 phase transition at t 苷 0. Snapshots (a) in an aggregation stage 共t 苷 10兲 and (b) at a later time 共t 苷 100兲, whereordering of the particles is observed. Each particle is accompanied by two 21兾2 charge point defects. Darkness presents the22value of qxx. Black and white correspond to qxx
苷 0 and 0.25,
respectively. Those correspond also to u 苷 0.25p, 0.75p, andu 苷 0, 0.5p, p as shown in the gradation map.

075502-4*Permanent address: Department of Physics, Kyoto University, Kyoto 606-8502 Japan.
Email address: ryoichi@scphys.kyoto-u.ac.jp[1] W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions (Cambridge University Press, Cambridge, 1995).
[2] V. M. Mostepanenko and N. N. Trunov, The Casimir Effectand its Application (Clarendon Press, Oxford, 1997).
[3] A. Borštnik, H. Stark, and S. Žumer, Phys. Rev. E 60, 4210
(1999).
[4] E. M. Terentjev, Phys. Rev. E 51, 1330 (1995).
[5] S. Ramaswamy, R. Nityananda, V. A. Raghunathan, andJ. Prost, Mol. Cryst. Liq. Cryst. 288, 175 (1996).
[6] R. W. Ruhwandl and E. M. Terentjev, Phys. Rev. E 56, 5561
(1997).
[7] H. Stark, Eur. Phys. J. B 10, 311 (1999); H. Stark, J. Stelzer,and R. Bernhard, Eur. Phys. J. B 10, 515 (1999).
[8] O. Mondain-Monval, J. C. Dedieu, T. Gulik-Krzywicki, andP. Poulin, Eur. Phys. J. B 12, 167 (1999).
[9] Y. G. Gu and N. L. Abbott, Phys. Rev. Lett. 85, 4719
(2000).
[10] P. Poulin, H. Stark, T. C. Lubensky, and D. A. Weitz, Science 275, 1770 (1997); P. Poulin, V. Cabuil, and D. A.
Weitz, Phys. Rev. Lett. 79, 4862 (1997); P. Poulin andD. A. Weitz, Phys. Rev. E 57, 626 (1998).
[11] R. W. Ruhwandl and E. M. Terentjev, Phys. Rev. E 55, 2958
(1997).
[12] T. C. Lubensky, D. Pettey, N. Currier, and H. Stark, Phys.
Rev. E 57, 610 (1998).
[13] B. I. Lev and P. M. Tomchuk, Phys. Rev. E 59, 591 (1999).
[14] V. A. Raghunathan, P. Richetti, and D. Roux, Langmuir
12, 3789 (1996); V. A. Raghunathan, P. Richetti, D. Roux,F. Nallet, and K. Sood, Langmuir 16, 4720 (2000).
[15] S. P. Meeker, W. C. K. Poon, J. Crain, and E. M. Terentjev,Phys. Rev. E 61, R6083 (2000); V. J. Anderson, E. M. Terentjev, S. P. Meeker, J. Crain, and W. C. K. Poon, Eur. Phys.
J. E 1, 11 (2001); V. J. Anderson and E. M. Terentjev, Eur.
Phys. J. E 1, 21 (2001).
[16] J. C. Loudet, P. Barios, and P. Poulin, Nature (London) 407,611 (2000).
[17] J. Yamamoto and H. Tanaka, Nature (London) 409, 321(2001).
[18] J. L. Billeter and R. A. Pelcovits, Phys. Rev. E 62, 711(2000).
[19] H. Löwen, P. A. Madden, and J. P. Hansen, Phys. Rev. Lett.
68, 1081 (1992); H. Löwen, J. P. Hansen, and P. A. Madden,J. Chem. Phys. 98, 3275 (1993).
[20] P. G. de Gennes and J. Prost, The Physics of Liquid Crystals (Clarendon, Oxford, 1993), 2nd ed.
[21] H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000).
[22] J. Fukuda and H. Yokoyama, Eur. Phys. J. E 4, 389 (2001).

075502-4