Document [original]

J. Chem. Phys. 135, 204702 (2011); https://doi.org/10.1063/1.3660377 135, 204702

Orientational prewetting of planar solid

substrates by a model liquid crystal

Cite as: J. Chem. Phys. 135, 204702 (2011); https://doi.org/10.1063/1.3660377

Submitted: 16 August 2011 • Accepted: 21 October 2011 • Published Online: 22 November 2011

Manuel Greschek and Martin Schoen

ARTICLES YOU MAY BE INTERESTED IN

Topological defects around a spherical nanoparticle in nematic liquid crystal: Coarse-grained

molecular dynamics simulations

The Journal of Chemical Physics 141, 114903 (2014); https://doi.org/10.1063/1.4894438

Simulation of Gay-Berne liquid crystal molecules confined to a plane

AIP Conference Proceedings 2220, 130016 (2020); https://doi.org/10.1063/5.0001132

Effects of flow on topological defects in a nematic liquid crystal near a colloid

The Journal of Chemical Physics 140, 054905 (2014); https://doi.org/10.1063/1.4862953

THE JOURNAL OF CHEMICAL PHYSICS 135, 204702 (2011)

Orientational prewetting of planar solid substrates by a model liquid crystal

Manuel Greschek1,a) and Martin Schoen1,2,b)

1Stranski-Laboratorium für Physikalische und Theoretische Chemie, Fakultät für Mathematik und

Naturwissenschaften, Technische Universität Berlin, Straße des 17. Juni 135, 10623 Berlin, Germany

2Department of Chemical and Biomolecular Engineering, Engineering Building I, Box 7905, North Carolina

State University, 911 Partners Way, Raleigh, North Carolina 27695, USA

(Received 16 August 2011; accepted 21 October 2011; published online 22 November 2011)

We present grand canonical ensemble Monte Carlo simulations of prewetting transitions in a model

liquid crystal at structureless solid substrates. Molecules of the liquid crystal interact via anisometric

Lennard-Jones potentials and can be anchored planar or homeotropically at the substrates. Fluid-

substrate attraction is modeled by a Yukawa potential of variable range. By monitoring the grand-

potential density and the nematic order parameter as functions of the chemical potential μ,sev-

eral discontinuous prewetting, wetting, and isotropic-nematic phase transitions are observed. These

transitions depend on both the range of the fluid-substrate attraction and the specific anchoring at

the substrate. Our results show that at substrates characterized by degenerate anchoring prewetting

occurs at lower μcompared with cases in which the anchoring is monostable. This indicates that

prewetting transitions are driven by orientational entropy because degenerate anchoring allows for

more orientationally distinct configurations of molecules compared with monostable anchoring. In

addition, by analyzing local density and various local order parameters, a detailed picture of the

[doi:10.1063/1.3660377]

I. INTRODUCTION

If a liquid crystal is in contact with a solid substrate, the

presence of that substrate may induce specific preferred orien-

tations of the molecules in its vicinity even though in the bulk

the liquid crystal would be isotropic under identical thermo-

dynamic conditions. The preferred orientation of molecules

in the immediate vicinity of the solid substrate is a conse-

quence of “anchoring”.1Anchoring may be perceived as an

energetic discrimination of a set of molecular orientations.

According to the definitions given by Jerôme,1anchoring con-

ditions may be classified as degenerate if an infinite number

of so-called “easy axes” exists, monostable if the molecules

align with a single preferred spatial direction, or multistable if

a finite number of orientations exists that are energetically dis-

criminated by a substrate. Perhaps, the earliest study of how

anchoring may affect properties of a liquid-crystal film has

been carried out by Sheng within the framework of Landau-

de Gennes theory.2,3Depending on the thermodynamic con-

ditions, Sheng observed an isotropic-nematic phase transition

occurring either near the solid substrate or far away from the

substrate.3Shortly after these early theoretical studies, it was

shown experimentally by Miyano4that liquid-crystal films

adsorbed on solid substrates may exhibit birefringence indi-

cating a significant degree of nematic order present in these

films. Understanding and eventually controlling the impact

of anchoring on wetting of solid surfaces is not only of aca-

demic interest but may also have important repercussions in

the realm of optoelectronic devises, for example.5

a)Electronic mail: [email protected].

b)Electronic mail: [email protected].

Therefore, it does not seem surprising that a host of ex-

perimental techniques exists to realize substrates with spe-

cific anchoring characteristics. According to Sonin, three

major groups can be identified: mechanical work, chemical

substance deposition, and external action.6For instance, a

substrate can be exposed to mechanical work by rubbing or

polishing it, by exposing it to photolithographic relief, or by

oblique evaporation of oxide films. Chemical substance depo-

sition can be realized either through deposition of surfactants

or polymeric films. The former causes homeotropic anchoring

whereas planar anchoring is realized by the latter approach.

Finally, external fields or flow can be used to manipulate the

anchoring of molecules at substrate surfaces as examples of

external action. Depending on details of how anchoring is

effected, molecules that are located at a larger/smaller dis-

tance from the substrate surface may be aligned in a specific

fashion. This, of course, requires a fluid-substrate interaction

potential of variable range. A more detailed discussion and

references to relevant literature can be found in Chapter 2 of

Sonin’s book.6

More recently, and as far as preparation of solid sub-

strates with specific anchoring is concerned, the focus has

shifted towards nanometer-size structures and heterogeneous

substrates. An example is the work by Yi et al. who use poly-

mer films nanoimprinted with checkerboard patterns of square

wells to demonstrate that the geometry of these wells allows

one to align liquid-crystal molecules in different ways.7They

confirm their experimental findings also by simple theoretical

considerations.

Once the substrate has been endowed with specific

anchoring characteristics the resulting structure of an ad-

jacent liquid crystal has most frequently been studied by

204702-2 M. Greschek and M. Schoen J. Chem. Phys. 135, 204702 (2011)

ellipsometry. For example, Hsiung et al.8have studied

the substrate-induced orientational order in an otherwise

isotropic liquid crystal and observe a weak nematic layer

on a silane-decorated glass substrate. The transition from

molecularly thin to thick films has been investigated by

Valignat et al. for droplets of nonvolatile liquids using

spatially resolved ellipsometry.9Chen et al. studied the ori-

entational wetting of a homologous series of cyanobiphenyl

liquid crystals on glass that was coated with silane surfactant

molecules.10 Lucht and Bahr used ellipsometry to investigate

the wetting behavior at the free surface of the isotropic

phase of binary mixtures of thermotropic liquid crystals.11

More recently, Bahr also reported a prewetting transition at

the interface between a nematic liquid crystal and water in

the presence of nonionic surfactant.12 In this work, it is also

demonstrated that increasing surfactant concentration drives

the prewetting transition towards a surface-critical point.

Experimental evidence for the existence of an orientational

prewetting transition is also provided through specific-heat

measurements carried out by Liu and co-workers.13

On the theoretical side, there has also been some interest

in wetting of solid substrates by liquid-crystalline materials.

Quite some time ago, Poniewierski and Sluckin have studied

wetting transitions.14 They employ a version of Maier-Saupe

mean-field theory properly modified for homeotropically an-

choring substrates. In their study, prewetting of the substrate

by a nematic phase is observed. Later, these same authors

used Landau-de Gennes’ theory to demonstrate that close to

the bulk isotropic-nematic phase transition the nematic film

near the substrate may be perceived as a low-temperature

Kosterlitz-Thouless phase.15 Landau-de Gennes theory has

been rather popular in the field of wetting phenomena.16–18

Another interesting example is the work by Patrício et al.19

These authors studied wetting transitions of nematic liquid

crystals on periodic, wedge-shaped substrates.

Using instead mean-field density functional theory Telo

da Gama investigated prewetting of solid substrates by liquid

crystals.20 If the fluid-substrate interaction decays exponen-

tially with increasing distance from the substrate, the prewet-

ting transition depends on how rapid this decay actually is.

Later, de las Heras and co-workers used a version of On-

sager’s theory21 to study phase transitions in a fluid com-

posed of hard spherocylinders confined to a slit-pore.22,23

In the limit of a large distance between the substrates in-

teresting wetting phenomena involving different mesophases

are observed. However, prewetting does not occur in this

system.

In addition, computer simulations have frequently been

employed to study orientational wetting of substrates by

liquid crystals. Dijkstra et al. investigated fluids composed

of hard spherocylinders interacting with hard substrates by

Monte Carlo (MC) simulation.24 By means of density func-

tional theory and MC simulation the same authors25 studied

wetting phenomena within Zwanzig’s26 simple hard-rod

model. MC simulations have been performed by Allen27

to investigate the impact of anchoring scenarios at solid

substrates on nematic liquid-crystalline phases. Again by

means of MC simulations, Downton and Allen investigated

the impact of various anchoring scenarios by grafting hard

spherocylinders homeotropically onto hard substrates.28 By

varying both the grafting density and the length of the grafted

molecules, a number of novel anchoring regimes in the

adjacent liquid-crystal phase could be observed. Using the

Gay-Berne model for a system of soft disks Bellier-Castella

et al. could demonstrate that the temperature of the nematic-

columnar transition depends on the substrate anchoring of the

disks.29 The Gay-Berne model has also been employed by

Müller et al. who observed that the formation of molecularly

thick films adsorbed on solid substrates does not involve a

prewetting transition under the conditions employed in their

simulations.30 By using a fluid composed of hard Gaussian

overlap particles interacting with planar substrates through

a hard-needle potential, Barmes and Cleaver were able to

realize both planar and homeotropic anchoring scenarios.31 In

their model, the substrate interacts with a “needle” of variable

length embedded inside a fluid molecule. If the length of the

needle corresponds to the long axis of a molecule, planar

alignment at the substrate results; if, on the other hand,

the needle shrinks to a point, homeotropic alignment is

observed.31 Cheung used MC simulations to demonstrate that

structured substrates may affect the alignment of molecules.32

Here, even tilted alignment of molecules could be realized by

fine-tuning the fluid-substrate interaction. Similarly, Zhang

et al. could show that a substrate may control the alignment of

molecules in (bulk-like) regions of an adjacent liquid-crystal

phase.33

In this paper, we will focus on the impact of various an-

choring scenarios on the prewetting transition and how this

transition depends on the range of the fluid-substrate inter-

action potential. We employ a simple model based upon an

anisotropic Lennard-Jones potential introduced by Hess and

Su34 and employed by us in previous studies.35–37 The remain-

der of this paper is organized as follows. In Sec. II,weintro-

duce the model on which our study is based. We summarize

the basic theoretical concepts in Sec. III. Results are presented

in Sec. IV and put into perspective in the concluding Sec. V.

Relevant general properties of the grand-potential density are

summarized in the Appendix.

II. MODEL SYSTEM

We consider a fluid composed of Nmolecules interacting

with each other in a pairwise additive fashion such that the

fluid-fluid (ff) configurational potential energy is given by

Uff(R,

U)=1



i=1



j=i

uff(rij ,

ui,

uj),(2.1)

where rij ≡ri−rjis the distance vector between the centers

of mass of particles iand jlocated at riand rj, respectively.

In Eq. (2.1),R≡{r1,r2,...,rN}and 

U≡{

u1,

u2,...,

uN}

are shorthand notations for the respective sets of center-of-

mass coordinates and unit vectors specifying the orienta-

tions of the liquid-crystalline molecules. Following earlier

work,34–39 we adopt the intermolecular interaction potential

204702-3 Orientational prewetting by a liquid crystal J. Chem. Phys. 135, 204702 (2011)

uff rij ,

ui,

uj=4σ

rij 12

−σ

rij 6

×{1+(

rij ,

ui,

uj)},(2.2)

where r=|r|and 

r=r/r. Throughout this work, notation

“...” indicates a unit vector. Hence, uff is a Lennard-Jones po-

tential where the attractive contribution is modified to account

for different relative orientations of a pair of molecules. In

Eq. (2.2),σdenotes the “diameter” of a spherical reference

molecule and is the depth of the attractive well in that ref-

erence model. The anisotropy of the fluid-fluid interaction is

accounted for by the function



rij ,

ui,

uj=5ε1P2(

ui·

uj)

+5ε2[P2(

ui·

rij )+P2(

uj·

rij )],

(2.3)

which one obtains from a summation of certain Wigner ma-

trices selected to preserve the head-tail symmetry of the

molecules34 and to guarantee that

d

uid

uj(

rij ,

ui,

uj)=0,(2.4)

which holds because in Eq. (2.3),

P2(x)=1

2(3x2−1),(2.5)

is the second Legendre polynomial. In other words, re-

mains unaltered if 

uiand/or 

ujare replaced by −

uiand/or

−

uj(head-tail symmetry).

In Eq. (2.3),ε1and ε2are (dimensionless) anisotropy

parameters determining the orientation dependence of inter-

molecular attractions. In this study, we deviate from our pre-

vious work35–37 and choose ε1=0.04 and ε2=−0.04. With

FIG. 1. Plots of u(r12,u1,u2) for a pair of molecules located in the x–y

plane such that r12 =(x12,y

12,0) and u1·u2=1. The white area at the

center of the plot is defined by the condition uff (r12,u1,u2)≥0 such that it

approximately represents the shape of the molecule at the center of the coor-

dinate system true to scale. In colored parts of the figure uff (r12,u1,u2)<0

where the color indicates the local value uff in units of (see attached color

bar).

this choice, molecules are nearly spherical as the plot in

Fig. 1reveals. Taking the ratio of zeros of uff in Eq. (2.2)

for side-side (

u1·

u2=±1, 

u1·

r=

u2·

r=0) and end-end

configurations (

u1·

u2=±1, 

u1·

r=

u2·

r=±1) as a def-

inition of the aspect ratio κof a fluid molecule, we obtain κ

≃1.095. The main purpose for selecting ε2=−0.04 is that

this reduces the overall fluid-fluid attraction compared with ε2

=−0.08 used in our previous works.35–37 As a consequence,

the density range of the bulk one-phase gas region is consid-

erably wider for temperatures that are not too much below the

(estimated) gas-(isotropic) liquid critical temperature.34

The fluid is confined between two planar, structureless

solid substrates separated by a distance sz0 along the z-axis

of a space-fixed Cartesian coordinate system. Hence, a fluid-

substrate (fs) contribution to the total configurational energy

maybecastas

Ufs(R,

U)=



k=1



i=1

u[k]

fs (ri,

ui),(2.6)

where

u[k]

fs (zi,

ui;η)=a1σ

zi±sz0/210

−a2

exp (−η|zi±sz0/2|)

|zi±sz0/2|g[k](

ui),

(2.7)

where zi=ri·

ezand 

ezpoints in the direction of the z-axis

of a space-fixed Cartesian coordinate system. We assume the

lower substrate (k=1) to be located at zw=−sz0/2 whereas

the upper one is located at zw=+sz0/2 (k=2). Because

the substrates are structureless, ufs at fixed molecular orienta-

tion depends only on the center-of-mass distance of molecule

ifrom either substrate along the z-axis. Moreover, because

we are interested in a physical situation in which a molecule

might interact with a single substrate at most we fix sz0 =50σ

throughout this work.

The orientation dependence of the fluid-substrate in-

teraction is described by the so-called anchoring function

0≤g[k](

ui)≤1 which we specify in Table I. We consider

combinations of anchoring functions that favor either planar

or homeotropic alignment. Notice that for nonspecific (n) an-

choring (g=1), no orientation of a molecule with respect

to the substrate plane is energetically favored. According to

the classification scheme of Jerôme1, this results in degener-

ate anchoring. Nevertheless, molecules prefer an on-average

homeotropic alignment. Notice that in our model, only the

center of mass of a molecule interacts with the substrate. This

is also the case in the model of Barmes and Cleaver if their in-

teraction “needle” embedded inside the elongated molecules

shrinks to a point (see Sec. I).31 In our model, the fluid-fluid

potential favors side-side arrangements of a pair of molecules

and stabilizes the homeotropic alignment at the substrate sur-

faces. In addition, because of κ>1 packing at the surface

plane is slightly more efficient if the weakly nonspherical

molecules are homeotropically rather than planar oriented. As

far as planar anchoring is concerned, the anchoring scenario

204702-4 M. Greschek and M. Schoen J. Chem. Phys. 135, 204702 (2011)

TABLE I. Combinations of anchoring functions g[k](u) used in this study

[see Eq. (2.7)]. Combinations of planar (p) and directional (d) as well as of

nonspecific (n) and homeotropic (h) anchoring are considered. Unit vectoreα

points along the α-direction of a space-fixed Cartesian coordinate system.

Acronym g[1](u)g[2](u)

pd (u·ex)2+(u·ey)2(u·ex)2

nh 1 (u·ez)2

referred to as directional (d) pins the molecules’ longer axes

to the x-direction whereas in the planar (p) scenario all ori-

entations on a unit circle in the substrate plane is energeti-

cally equivalent. Again, according to Jerôme’s classification

scheme, the former causes monostable anchoring whereas in

the latter case anchoring is degenerate.

In Eq. (2.7), we employ a Yukawa-like attractive contri-

bution, where η−1may be interpreted as a “screening length”

which determines the range of the fluid-substrate attraction.

Notice that the Yukawa-type of attraction—usually describ-

ing the interaction between screened charges—has been cho-

sen without any specific experimental system in mind. Rather

this form of fluid-substrate attraction has been selected out

of mere convenience because it permits us to tune the range

of this attraction continuously by varying η. In this work, we

consider ησ =0.25, 0.45, and 1.0 such that with increasing

ηthe fluid-substrate attraction becomes increasingly shorter

range. Dimensionless parameters,

a1≡−

ησ +1

ησ −9,(2.8a)

a2≡−

exp (ησ)

ησ −9,(2.8b)

have been introduced to guarantee that

du[k]

dzzmin

=0,(2.9a)

u[k]

fs (zmin)=−

4, (2.9b)

irrespective of ησ, where zmin =±sz0/2 ∓σdenotes the lo-

cation of the minimum of ufs relative to the upper and lower

substrate surface, respectively. The depth of the attractive well

at z=zmin in Eq. (2.2) has been chosen arbitrarily but large

enough to guarantee complete wetting of the substrates as

one approaches gas-(isotropic) liquid coexistence. We illus-

trate the various fluid-substrate potentials in Fig. 2. The plots

show that with our choice of a1and a2only the range of

the attractive part of the total fluid-substrate attraction is af-

fected by the specific choice of ησ. From the inset in Fig. 2,

it is also clear that molecules located near the midpoint z=0

are nearly unaffected by the presence of either wall because

u[1]

fs +u[2]

fs ≈0 irrespective of the value of ησ considered in

this study. In other words, portions of the liquid crystal in the

vicinity of one wall are indeed not subject to a simultaneous

interaction with the other one.

-4.00

-3.00

-2.00

-1.00

0.00

-0.50 -0.45 -0.40 -0.35 -0.30

ufs

z/sz0

-0.06

-0.04

-0.02

0.00

-0.30 -0.20 -0.10 0.00

FIG. 2. Total fluid-substrate interaction potential u[1]

fs +u[2]

fs as a function of

position zin the lower half of the system and for η=0.25 (– ·–·), η=0.45

(—), and η=1.00 (----) [see Eq. (2.7)]. Inset shows an enhancement close to

the midpoint of the system. Data for η=1.00 are omitted because they are

indistinguishable from u[1]

fs +u[2]

fs =0.

III. THEORETICAL BACKGROUND

A. Key concepts of phenomenological and statistical

thermodynamics

In this paper, we are concerned with formation and ther-

modynamic stability of liquid-crystalline films adsorbed on

solid surfaces. Therefore, it seems most appropriate to con-

sider the grand canonical ensemble which permits relatively

easy access to the grand potential . In the grand canonical

ensemble

(μ, T, A, sz)=−kBTln (μ, T, A, sz),(3.1)

where μdenotes the chemical potential, Tis the temperature,

kBis the Boltzmann’s constant,

=

NIexp(βμ)

m5N

Z,(3.2)

is the partition function of the grand canonical ensemble

(β≡1/kBT), and

Z=1

2NN!dRd

Uexp(−βU),(3.3)

is the configuration integral.35 In Eq. (3.2),≡βh2/2πm

is the thermal de Broglie wavelength, mis the molecular mass,

and Iis the moment of inertia of a molecule. The prefactors

1/2Nand N!inEq.(3.3) correct for double counting equiv-

alent configurations that arise by replacing the set 

Uby the

equivalent set −

U(see Sec. II) and take notice of the in-

distinguishability of the molecules, respectively. The quantity

U=Uff +Ufs in the integrand is the total configurational po-

tential energy of the liquid crystal [see Eqs. (2.1) and (2.7)].

By purely thermodynamic reasoning,40,41 one can show

that35

d=−SdT−Ndμ−Psz0dA−PzzA0dsz,(3.4)

where Sdenotes entropy, Nis the number of molecules, and

P≡1

2(Pxx +Pyy) and Pzz are the respective transverse and

normal (diagonal) components of the pressure tensor P(with

204702-5 Orientational prewetting by a liquid crystal J. Chem. Phys. 135, 204702 (2011)

respect to the position of the substrate plane). In the absence

of shear strains Pcan be represented by the 3 ×3matrix

P=⎛

⎜

⎝

P00

0P0

00Pzz

⎞

⎟

⎠.(3.5)

In writing Eq. (3.4), we assume components of Pto act on the

faces of a rectangular fluid lamella where A0=sx0sy0 such

that sα0is the side length of the lamella in the direction of the

α-axis of a Cartesian coordinate system.40 Subscript “0” indi-

cates that these side lengths refer to dimensions of the lamella

in some unstrained reference state. Consequently, Aand szin

Eq. (3.4) correspond to substrate area and separation but in a

strained (i.e., compressed or dilated) state. This state can be

realized in a gedankenexperiment by moving the boundaries

separating the lamella from its surroundings.40,41 Hence, the

last two terms in Eq. (3.4) account for mechanical work ex-

changed between the lamella and the surroundings.

The advantage of this approach is that maybeper-

ceived as a homogeneous function of degree one in Abecause

ufs depends only on the distance between a molecular center-

of-mass and the substrate plane along the z-axis. Hence, ufs

is translationally invariant in both x- and y-directions.40 One

may therefore apply Euler’s theorem and obtain

ω≡

V=−P,(3.6)

for the grand-potential density, where V=Asz0 is the vol-

ume of the lamella. Employing the so-called hypervirial

theorem,42 it follows from Eqs. (2.2) to (2.5) that35

Pαα =ρkBT−24

V

(i,,j)σ

rij 12

−1

2σ

rij 6

×(

rij ,

ui,

uj)rα

ij rα

ij 

−30εε2

V

(i,j)σ

rij 6

ui·

rij uα

irα

+

uj·

rij uα

jrα

ij −(

ui·

rij )2+(

uj·

rij )2

×rα

ij rα

ij ,α=x,y. (3.7)

On the right side of Eq. (3.7), superscripts αrefer to the α-

component of the associated vector and angular brackets de-

note a grand canonical ensemble average. The first term of

the right side is the ideal-gas contribution where ρ≡N/V

is the number density. Together with Eq. (3.6),Eq.(3.7) per-

mits a quantitative discussion of phase equilibria. Additional

properties of the grand-potential density are discussed in the

Appendix.

B. Properties

Besides the grand-potential density already introduced in

Eq. (3.6) measures of nematic order and structure of the ad-

sorbed liquid crystals will be considered in this section. To

that end, we introduce the so-called (instantaneous) alignment

tensor43

Q=1



i=1

(3

ui⊗

ui−1),(3.8)

which is a second-rank tensor that can be represented by a

3×3 matrix. In Eq. (3.8),1denotes the unit tensor and

the operator “⊗” stands for the direct (i.e., “dyadic”) prod-

uct. From Eq. (3.8), it is also clear that the alignment ten-

sor is real and symmetric. Moreover, it is traceless such that

Tr Q=0 where the operator “Tr ” represents the trace. Using

Jacobi’s method,44 one may compute the eigenvalues λ−≤λ0

≤λ+of Qand the associated eigenvectors. Following stan-

dard practice,35,36,45–49 we define the nematic (Maier-Saupe)

order parameter as S≡λ+. The eigenvector associated with

the instantaneous eigenvalue λ+is taken as the nematic di-

rector 

n. Moreover, Qmay be diagonalized in the basis of its

three eigenvectors to obtain diag Q. A basic theorem of linear

algebra establishes the trace of a tensor as one of its scalar

invariants50 such that Tr Q=Tr (diag Q)=0,where

diag Q=⎛

⎜

⎝

−λ+/2−ζ00

0−λ+/2+ζ0

00λ+

⎞

⎟

⎠,(3.9)

which was demonstrated by Low who uses an expansion in

terms of Wigner matrices.51 In Eq. (3.9),ζis a measure of

apparent biaxiality. As Eppenga and Frenkel already pointed

out, ζis nonzero in any finite system even in the isotropic

phase.47 As we have recently demonstrated the associated os-

tensible biaxiality order parameter ξ≡ζmaybeusefulin

locating the isotropic-nematic (IN) phase transition in a rea-

sonably large system.37

In addition to these global measures of order, it will also

be interesting to determine order locally. To that end one no-

tices that on account of the fluid-substrate interaction [see

Eq. (2.7)], 

uiin Eq. (3.8) may be perceived as a local vari-

able 

ui(z). This, of course, permits to define a local alignment

tensor Q(z) such that a local nematic order parameter S(z)

=λ+(z)may be introduced. However, some caution is ad-

visable at this point. In computer simulations, S(z) may be

computed as a histogram by considering the orientations of

N(z) molecules located within a small volume Aδzcentered on

z. If this volume is relatively small so will be N(z). Because

of Eq. (2.2), these molecules tend to align their longer axes

with some local director such that there is always some artifi-

cial local nematic order even though the entire liquid crystal

is still in the isotropic phase.47 To correct for this finite-size

effect Richter and Gruhn45 devised an efficient scheme that

we employ throughout this work.

In addition, another measure of local order will turn out

to be useful. Its definition is based upon the realization that

the z-axis is a distinguished direction in our system due to the

presence of the planar substrates. Hence, we define

P2(z)=1

21

N(z)



i=13[

ui(z)·

ez]2−1.(3.10)

204702-6 M. Greschek and M. Schoen J. Chem. Phys. 135, 204702 (2011)

From this definition, it is immediately clear that P2(z)→−

if molecules are oriented in a parallel fashion with respect to

the substrate plane; P2(z)→1, if they are homeotropically

aligned instead. In this latter case, S(z)=P2(z) whereas S(z)

and P2(z) will generally deviate as P2(z)→−

2. Notice also

that S(z)≥0, where the equality holds if and only if Q(z)=0,

where 0denotes the zero tensor. In this latter case, λ−=λ0

=λ+=0 indicating an ideal isotropic arrangement of the

molecules.

Finally, we introduce the local density of the fluid ρ(z)

as a measure of probability to find the center of mass of a

molecule at a position z. Because of the form of the fluid-

substrate potential introduced in Eq. (2.7), the local density

must be translationally invariant with respect to the x- and y-

components of the center of mass position r. Formally, one

may introduce ρ(z) via the expression

ρ(z)≡N



i=1

δ(z−zi)=N(z)

Aδz ,(3.11)

where δ(z−zi)istheDiracδ-function and the far right side

of the equation gives an operational expression useful for the

computation of ρ(z) in a computer simulation.

IV. RESULTS

A. Numerical details

In this work, we employ MC simulations in the grand

canonical ensemble employing the algorithm discussed in de-

tail by Gruhn and Schoen.52 This algorithm closely resembles

that proposed by Adams for a simple Lennard-Jones fluid.53

It allows one to realize numerically a Markov process that

generates a distribution in configuration space proportional to

exp{−β[U(R,

U)−μN]−ln N!−5Nln(m/I)}in agree-

ment with Eqs. (3.2) and (3.3).

We refer to a MC cycle as a sequence of Nattempts to

displace or rotate a molecule and Nattempted creations of

new or removals of already existing molecules where Nis

the actual number of molecules present at the beginning of a

new cycle. To avoid biasing, the generation of configurations,

displacements, and rotations as well as creation and removal

are attempted with equal probability. Our simulations are

based upon 2 ×104cycles for equilibration followed by 2

×105cycles to compute ensemble averages. As we have

verified for a couple of thermodynamic states, these numbers

are large enough to guarantee reliable results even in the

immediate vicinity of a phase transition. To save computer

time, we employ a combination of a conventional Verlet

with a link-cell neighborlist as described in Allen and

Tildesley’s book.54 In addition, fluid-fluid interactions are cut

off beyond an intermolecular separation of rc=3σ, which

we use throughout this work; no such cutoff is applied to

fluid-substrate interactions.

We begin a sequence of simulations at a sufficiently low

μat which a bulk system would be a thermodynamically sta-

ble, dilute gas. Subsequent runs are started from the final con-

figuration of the preceding one where we increase μin small

steps δμ. During this sequence, we monitor ωcomputed via

Eqs. (3.6) and (3.7). At some characteristic value of μ,the

system will become thermodynamically unstable and under-

goes a transition from phase ito some other phase j.How-

ever, this transition does not necessarily occur at coexistence,

that is the associated grand-potential densities do not neces-

sarily satisfy Eq. (A5) but the inequality ωj<ω

iinstead. In

other words, the transition typically takes the system from a

metastable to a thermodynamically stable phase. To locate μij

[see Eq. (A5)], one can start from the current state point in

phase jand decrease μwhich causes the associated ωjto in-

crease because of Eq. (A1) until Eq. (A5) is satisfied.

This procedure is valid for phase transitions involving a

sufficiently large density difference between phases iand j

(for example, a transition from a gas-like film to a liquid-like

phase). However, for reasons to be explained below, we are

also interested in locating the IN phase transition which for

the present model is accompanied by a rather minute density

change.35,39 In this case, it will be useful to monitor the vari-

ation of the nematic order parameter Swith μ.

Throughout this work, we express all quantities of inter-

est in terms of the customary dimensionless (i.e., “reduced”)

units. For example, length is given in units of σ, energy in

units of , and temperature in units of /kB. Other derived

quantities are expressed in terms of suitable combinations of

these basic quantities. For example, pressure-tensor elements

are given in units of /σ3. In all the simulations, we consider

a computational cell of sx0 =sy0 =20, sz0 =50 to guarantee

that our model system effectively mimics a fluid phase inter-

acting with a single substrate as we already demonstrated in

Sec. II. Histograms of ρ(z), S(z), and P2(z) are computed with

abinwidthδz=0.1. For this system size and over the range

of chemical potentials studied between 900 and 17 000 parti-

cles are accommodated on average at a constant temperature

T=0.9 below the bulk gas-liquid critical point.34

B. Phase transitions

We begin our discussion of results with plots of ω

versus μin Fig. 3for nh anchoring. For the shortest-range

fluid-substrate potential characterized by η=1.00, the plot

in Fig. 3(a) shows that a low-density phase is thermodynam-

ically stable for μ−13.00. At this μ, the plot of ω(μ)for

the low-density phase intersects another curve ω(μ) with a

slope significantly larger in magnitude corresponding to a

higher-density phase [see Eq. (A1)]. The precise physical

nature of these different phases requires an analysis of ρ(z),

S(z), and P2(z), which we defer to Sec. IV C. Beyond the

intersection (i.e., as μincreases further), the low-density

phase remains metastable over a certain range of chemical

potentials and apparently becomes unstable at that chemical

potential at which the corresponding ω(μ) curve ends. The

same statements apply to ω(μ) for the high-density phase

below the chemical potential of intersection. If one increases

the range of the fluid-substrate interaction potential to η

=0.45, plots in Fig. 3(b) show that a third phase of inter-

mediate density arises such that two chemical potentials

exist at which pairs of phases coexist. Increasing the range

of fluid-substrate interactions by decreasing the screening

204702-7 Orientational prewetting by a liquid crystal J. Chem. Phys. 135, 204702 (2011)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(a)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(a)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(a)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(b)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(b)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(b)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(b)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(c)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(c)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(c)

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

(c)

FIG. 3. Plots of grand-potential density ωand (global) nematic order param-

eter Sas functions of chemical potential μfor nh anchoring (see Table I).

(a) η=1.00 [see Eq. (2.7)], (◦)ω,()Sfor a gas-like adsorbed film and

(•)ω,()Sfor a liquid-like phase. (b) As (a) but for η=0.45; notice that

a prewetting transition occurs where (+)representsωand (×) refers to S.

dynamic stability for specific phases whereas solid lines are fits to the discrete

data points to guide the eye.

constant to η=0.25, the plot in Fig. 3(c) reveals that the

one-phase region for the intermediate-density phase shifts to

slightly lower chemical potentials and widens somewhat.

The corresponding plot of the global nematic order pa-

rameter Sin Fig. 3(a) indicates that with increasing chemical

potential nematic order increases slightly in the low-density

phase but remains relatively small even at phase coexistence.

Hence, the low-density phase remains more or less isotropic

globally. At phase coexistence, Sdrops slightly and discontin-

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

-13.50 -13.25 -13.00 -12.75 -12.50 -12.25

0.0

0.2

0.4

0.6

0.8

1.0

FIG. 4. As Fig. 3(c), but for the pd anchoring scenario (see Table I). Inset is

an enlargement of the part of the plot within the box.

uously indicating that the newly formed high-density phase is

more isotropic. Interestingly, when the intermediate-density

phase forms in Figs. 3(b) and 3(c), nematic order increases

substantially reflected by an increase of Sby roughly a factor

of 2. When the intermediate-density phase gives way to the

high-density phase, this high nematic order is lost at first.

However, in all parts of Fig. 3, one realizes that eventually

there is a strong increase in nematic order at sufficiently high

μ≃−12.50. This value of μis almost independent of the

range of the fluid-substrate interaction potential indicating

that nematization is mostly occurring in the bulk of the

high-density phase, that is in regions that are not subject to

strong fluid-substrate interactions. However, the shape of the

plot of Sversus μbecomes a bit more rounded as ηdecreases.

The formation of the nematic phase is accompanied by

a small change in mean density of the liquid crystal.35 Be-

cause of Eqs. (A1) and (A3), one would therefore anticipate

a change both in slope and curvature in the plots of ω(μ)at

that μat which the IN phase transition occurs. However, this

effect seems to be much smaller than the statistical accuracy

with which the curves ω(μ)inFig.3can be obtained.

These general features are retained if the pd instead of

the nh anchoring scenario is considered. Therefore, to save

space we show results obtained for the pd anchoring only for

the case η=0.25 in Fig. 4as a typical illustration. The most

significant difference between the plots in Fig. 4and the cor-

responding Fig. 3(c) concerns the existence of four instead of

three phase transitions over the range of chemical potentials

considered. Increasing μfrom the smallest value for which

data have been plotted in Fig. 4, a low-density phase of nearly

constant ωis thermodynamically stable until μ≃−13.28.

Over this range of μ’s, the corresponding plot of Sshows

that there is relatively little nematic order in our system. At

μ≃−13.28, two of the curves ω(μ) of slightly different slope

intersect such that a phase of slightly larger mean density is

thermodynamically stable for μ−13.28. This new phase

exhibits a somewhat enhanced degree of nematic order but S

<0.2 still remains relatively small. At μ≃−13.18 a second

discontinuous phase transition arises. From the inset in

Fig. 4, one notices that in the newly formed phase the mag-

nitude of the slope of ω(μ) is again larger. Hence, this new

204702-8M. Greschek and M. Schoen J. Chem. Phys. 135, 204702 (2011)

FIG. 5. “Snapshots” from MC simulations in the x–zplane. The solid sub-

strates (not shown) are at the left and right ends of the rectangular slabs, re-

spectively. The black solid line represents the simulation box. The snapshots

correspond to η=0.25 and μ=−12.25, where the system in both nh (upper

plot) and pd anchoring scenarios (lower plot) is deep in the nematic phase

[see Figs. 3(c) and 4]. Choosing the horizontal direction (z-axis) particles are

colored in blue if ui·ez=1,whereas particles in red satisfy ui·ez=0.

phase is characterized by an elevated mean density compared

with its immediate predecessor at lower μ. An inspection

of the corresponding data for Sshows that nematic order

increases quite substantially at the second phase transition but

remains smaller than S=0.4 over the range of μ’s for which

this phase remains thermodynamically stable. This range

ends at μ≃−13.02 such that for μ−13.02 a high-density

phase characterized by the largest magnitude of the slope of

ω(μ) is thermodynamically stable. When this phase forms at

μ≃−13.02, Sdrops by almost a factor of 2 such that this

high-density phase is mostly isotropic.

Eventually, at μ≃−12.50 an IN phase transition occurs

and Srises steeply as μincreases further. A comparison with

plots in Fig. 3(c) shows that the IN phase transition again

occurs at about the same μirrespective of the anchoring

scenario. However, the direction of 

nis determined by the

specific anchoring as plots in Fig. 5illustrate. The upper

plot in Fig. 5is a representative “snapshot” from the MC

simulation illustrating the structure of the nematic phase for

the nh anchoring scenario. As one can see molecules are

mostly homeotropically aligned with a considerable number

of defects (i.e., particles with orientations other than those

imposed by the specific anchoring scenario). Because the

fluid-substrate potential at the center of the simulation box

has approximately decayed to zero (see Fig. 2) the dominant

homeotropic alignment is mediated by molecules that are

closer to the substrate. Hence, it is possible to “imprint” the

chemical nature of the substrates represented by a specific

anchoring scenario onto an adjacent liquid-crystalline phase

through a collective effect and over relatively large distances.

A more quantitative discussion of local structural features of

nematic phases and different anchoring scenarios is, however,

postponed until Sec. IV C.

One also notices from Figs. 3(c) and 4that the scatter of

Sin the high-density isotropic phase (−13.00 ≤μ≤−12.5)

is substantially larger for the pd compared with the nh anchor-

ing scenario which we rationalize as follows. At the lower

substrate the anchoring function g[1](

u) is degenerate in the

sense of Jerôme1because molecules do not suffer any energy

penalty as long as they align their long axes with any direction

on the unit circle in the x–yplane. On the other hand, g[2](

causes monostable (d) anchoring. If Sis relatively small, only

molecules in the immediate vicinity of the substrate exhibit

a certain degree of local nematic order with associated local

directors 

n[1](z) and 

n[2](z). These local nematic regions are

separated by a relatively large isotropic regime. Because

of its weakness, local nematic order is short range. Due to

monostable anchoring, 

n[2](z) remains more or less constant.

On the contrary, 

n[1](z) may vary in the x–yplane because of

the degenerate planar anchoring. Whereas this does not affect

the local degree of nematic order, any angle between 

n[1](z)

and 

n[2](z) affects the global order parameter computed via

Qwhich involves a sum over all molecular orientations [see

Eq. (3.8)]. This argument also implies that once the system

is in the nematic phase, molecules local and global directors

become fixed in space such that the scatter in Sshould be

smaller. Indeed, our results for μ−12.50 displayed in

Fig. 4comport with this line of arguments. Moreover, we

realize that none of the plots of Sin Fig. 3suffers from similar

statistical accuracies in the isotropic high density regime. This

is because at the substrate molecules align always homeotrop-

ically regardless of whether they are anchored nonspecifically

or homeotropically.35,36 In these cases, similar reorientation

effects are therefore prevented by both anchoring functions.

C. Structure and local orientation

The question now becomes: what is the physical nature

of the different phases represented by the plots in Fig. 3?

This question may be addressed directly by considering lo-

cal structure and orientation of the various phases. We be-

gin discussing these features in Fig. 6for the shortest-range

fluid-substrate potential characterized by η=1.00. As the

corresponding plots in Fig. 3(a) indicate two phases coexist at

μ≃−13.00. Corresponding plots of ρ(z)inFig.6(a) show

that the lower-density phase consists essentially of a thin film

of three fluid layers at the nonspecifically anchoring, lower

substrate (z/sz0 =−0.5) whereas only a monolayer of very

low density is adsorbed at the homeotropically anchoring up-

per substrate (z/sz0 =+0.5). These layers are reflected by pe-

riodic maxima of ρ(z). As we have verified, lowering μcauses

molecules to desorb continuously from the nonspecifically an-

choring substrate until only a submonolayer film is left.

Corresponding plots of both S(z) and P2(z) reveal that

molecules are homeotropically aligned where the degree of

homeotropic alignment is larger at the upper substrate as

one would have anticipated. At the nonspecifically anchor-

ing, lower substrate one sees that S(z)≈P2(z) as far as the

first three layers are concerned. However, as one approaches

204702-9 Orientational prewetting by a liquid crystal J. Chem. Phys. 135, 204702 (2011)

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(a)

0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

S(z),P2(z)

z/sz0

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(b)

0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

S(z),P2(z)

z/sz0

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(c)

0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

S(z),P2(z)

z/sz0

FIG. 6. Plots of local density ρ(z) (—) and local orientational order param-

eters P2(z)(—)andS(z) (—) as functions of position zbetween lower (z/sz0

=−0.5) and upper substrate (z/sz0 =+0.5) (sz0 =50). Data are obtained

for nh anchoring and η=1.00. Plots in parts (a) and (b) have been obtained

for the same chemical potential μ=−13.00 and correspond to coexisting

phases of lower and higher mean density, respectively. The plots in part (c)

have been obtained at μ=−12.20 deep in the nematic high-density phase

[see Fig. 3(a)].

the interface between the thin film and the adjacent bulk-like

gas phase S(z)= P2(z) indicating that molecules prefer to

align themselves parallel with the substrate plane. The pre-

ferred alignment at the film-gas interface is signaled by a pro-

nounced peak in S(z) at about z/sz0≃−0.42 in Fig. 6(a).Note

also that this peak in S(z) is not an artifact caused by the rela-

tively small number of molecules located at the film-gas inter-

face because our data have been corrected for such finite-size

effects as we explained in Sec. III B.

The structure of the coexisting phase can be analyzed

through the plots displayed in Fig. 6(b). Comparing these lat-

ter plots with the ones in Fig. 6(a), one notices that the liquid

crystal has a bulk-like, homogeneous region where ρ(z)=ρb

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(a)

0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

S(z),P

2(z)

z/sz0

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(b)

0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

S(z),P2(z)

z/sz0

FIG. 7. As Fig. 6,butforη=0.25. Plots in parts (a) and (b) have been

obtained for the same chemical potential μ=−13.30 and correspond to

coexisting phases of lower and intermediate mean density, respectively [see

Fig. 3(c)].

and ρbis the density of a corresponding bulk liquid crystal

at the same Tand μ. Closer to each substrate the molecules

arrange their centers of mass in individual layers parallel with

the substrate plane. At the nonspecifically anchoring sub-

strate, nematic order remains relatively low compared with

the homeotropically anchoring substrate. As one departs from

the substrate planes, nematic order decays such that the bulk-

like region of the liquid crystal remains isotropic. In the

more ordered portions close to the substrates, molecules are

homeotropically aligned as revealed by the relation P2(z)

≃S(z). Notice also that in Fig. 6(b), the specific anchoring

scenario affects only the orientation of the molecules but not

their packing. This is inferred from the plot of ρ(z)inFig.6(b)

which, unlike its counterpart in Fig. 6(a), remains symmetric

with respect to the midplane located at z=0.

If the fluid-solid interaction potential becomes longer

range plots in Figs. 3(b) and 3(c) indicated that an

intermediate-density phase arises which is accompanied by

a sudden increase in global nematic order. Plots in Figs. 7

and 8illustrate the structure of the coexisting phases at μ=

−13.30 and −13.05 and for η=0.25. As shown before in

Fig. 6(a), the corresponding plot in Fig. 7(a) shows that a

film of three distinct layers is adsorbed at the nonspecifically

anchoring substrate whereas a very dilute monolayer is sta-

bilized by the homeotropically anchoring substrate. Orienta-

tional effects in Fig. 7(a) are also very similar to those visible

from the corresponding plots in Fig. 6(a). This applies in par-

ticular to the maximum in S(z) at the film-gas interface located

approximately at z/sz0≃−0.4. However, the corresponding

204702-10 M. Greschek and M. Schoen J. Chem. Phys. 135, 204702 (2011)

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(a)

0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

S(z),P2(z)

z/sz0

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(b)

0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

S(z),P

2(z)

z/sz0

FIG. 8. As Fig. 6,butforη=0.25. Plots in parts (a) and (b) have been

obtained for the same chemical potential μ=−13.05 and correspond to co-

existing phases of higher and intermediate mean density, respectively [see

Fig. 3(c)].

plot in Fig. 7(b) reveals that at the homeotropically anchor-

ing substrate the adsorbed film has undergone a prewetting

transition from essentially a thin submonolayer to a micro-

scopically thicker film of three distinct layers. This thicker

film also exhibits distinct nematic order because the condi-

tion P2(z)=S(z) is satisfied until the film-(bulk) gas inter-

face is reached. Its degree of nematic order is substantially

larger than that observed for the film adsorbed at the opposite,

nonspecifically anchoring substrate as one would have antici-

pated. However, these directions are slightly different for the

film adsorbed at the nonspecifically and that at the homeotrop-

ically anchoring substrate. As one can see from the corre-

sponding plots of P2(z), this quantity passes through a weak

minimum when S(z) becomes maximum at the interface at the

nonspecifically anchoring substrate. At the homeotropically

anchoring substrate, both S(z) and P2(z) remain positive at the

interface. Similar features have been reported earlier at the

nematic-gas interface of homogeneous, free-standing liquid

crystal films.55,56 Hence, as far as nh anchoring is concerned

an already ordered submonolayer film undergoes a prewetting

transition to form a molecularly thick nematic film.

As the chemical potential increases over the range

−13.30 ≤μ≤−13.05, the films adsorbed at both substrates

increase in thickness. This can be seen by comparing plots in

Figs. 7(b) and 8(b) which show that for both films the film-gas

interface has shifted to smaller values of |z|. Notice that at

both substrates, one additional layer of molecules is visible in

the plot of ρ(z)inFig.8(b) compared with the corresponding

one in Fig. 7(b).Atμ=−13.05, the gas phase separating the

adsorbed films at both substrates condenses spontaneously

such that again a bulk-like, homogeneous fluid forms as

the plot in Fig. 8(b) shows. Comparing Fig. 8(b) with the

corresponding plot in Fig. 6(b), one notices that nematic order

induced by both substrates is much more pronounced if the

fluid-substrate interaction is longer range. Moreover, as one

approaches the midplane at z=0, the preferential orientation

of the molecules extends much farther into the bulk-like

region of the liquid crystal. One also notices that in this bulk-

like region ρ(z) does not correspond to ρbbut decays slightly

towards the midplane. This is because the slow decay of

orientational order as one departs from either substrate causes

the packing of molecules to gradually change. This effect

has been analyzed earlier by us36 for the same model liquid

crystal but for molecules with a slightly larger aspect ratio.

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(a)

0.3 0.4 0.5

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

S(z),P

2(z)

z/sz0

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(b)

0.3 0.4 0.5

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

S(z),P2(z)

z/sz0

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(c)

0.3 0.4 0.5

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

S(z),P

2(z)

z/sz0

FIG. 9. As Fig. 7, but for the pd anchoring scenario. (a) μ=−13.29, (b) μ

=−13.28, and (c) μ=−13.18 (see also Fig. 4).

204702-11 Orientational prewetting by a liquid crystal J. Chem. Phys. 135, 204702 (2011)

Finally, one notices from Fig. 3(c) that at sufficiently

large μ,theglobal order parameter for the high-density,

liquid-like system increases. This indicates again that the liq-

uid crystal undergoes an IN phase transition. The structure of

the nematic phase at η=0.25 turns out to be very similar to

what has already been illustrated by the plots in Fig. 6(c) for

the shorter-range fluid substrate potential (η=1.00). In the

interest of conciseness, we therefore refrain from presenting

corresponding data for η=0.25.

At this stage, it seems interesting to investigate the struc-

ture of coexisting phases for the pd anchoring scenario.

We are again restricting this discussion to the longest-range

fluid substrate potential (η=0.25). At sufficiently low μ

−13.28, thin submonolayer films are adsorbed at both sub-

strates [see Fig. 9(a)]. Here, molecules have a preferred paral-

lel orientation with respect to the substrate planes as expected.

This is inferred from a relatively large value of S(z) whereas

P2(z) remains negative. The submonolayer adsorbed at the

lower, parallel anchoring substrate is in equilibrium with a

thicker film comprising roughly four layers of molecules.

Hence, at μ≃−13.28, a prewetting transition occurs at the

planar anchoring substrate [see Fig. 9(b)]. The submonolayer

film adsorbed at the directionally anchoring upper substrate

remains thermodynamically stable until at μ≃−13.18, a

second prewetting transition occurs during which the mono-

layer at the directionally anchoring substrate is transformed

into a thicker film of about four distinct layers as the plots

in Fig. 9(c) illustrate. Again, molecules pertaining to these

films exhibit a preferred planar orientation with respect to the

substrate plane. As μincreases, the films at both substrates

grow in thickness continuously until at μ≃−13.02, the entire

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(a)

0.3 0.4 0.5

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

S(z),P2(z)

z/sz0

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

(b)

0.3 0.4 0.5

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

S(z),P2(z)

z/sz0

FIG. 10. As Fig. 9,butforμ=−13.02 (a) and μ=−13.00 (b).

0.0

0.5

1.0

1.5

2.0

2.5

-0.5 -0.4 -0.3

ρ(z)

0.3 0.4 0.5

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

S(z),P

2(z)

z/sz0

FIG. 11. As Fig. 9,butforμ=−12.25 (see also Fig. 4).

system condenses and a mostly isotropic, high-density phase

forms (see Fig. 4).

A comparison of plots of S(z)inFigs.9(c) and 10(a) show

that at the lower, planar aligning substrate, the order of the

adsorbed film increases only slightly with μ. At the upper,

directionally aligning substrate, the degree of local order re-

mains nearly the same as the film thickens. Local nematic or-

der at the lower substrate increases only in the presence of a

stabilizing isotropic liquid as a comparison between plots in

Figs. 10(a) and 10(b) reveals. This observation is in contrast

with findings by Rull et al. who observed surface nematiza-

tion in a thin film without the presence of an adjacent isotropic

liquid.57 However, we would see an effect similar to the one

reported by Rull et al. but only for already metastable nematic

films.

Increasing μfurther causes an IN phase transition in the

high-density phase at μ≃−12.50 as the corresponding plots

in Fig. 4clearly show. The structure of this high-density ne-

matic phase is illustrated by the plots in Fig. 11. As one can

see there molecules in the entire system have a preferred pla-

nar orientation with the solid substrates because S(z)islarge

whereas P2(z) is negative. The structure of the nematic phase

is also illustrated by the “snapshot” from an MC simulation

presented as the lower part of Fig. 5.

V. DISCUSSION AND CONCLUSIONS

In this work, we employ MC simulations in the grand

canonical ensemble to study wetting of solid substrates by a

model liquid crystal. We focus on different anchoring scenar-

ios in which specific orientations of the molecules are ener-

getically discriminated. We consider two types of anchoring

scenarios that can be described as “degenerate” and “mono-

stable”. Anchoring scenarios pertaining to the latter group

(i.e., homeotropic and directional anchoring) favor orienta-

tions of molecules along a single easy axis whereas those

scenarios belonging to the former class give molecules the

freedom to orient themselves either in all three spatial direc-

tions (i.e., nonspecific anchoring) or in any direction within

the substrate plane (i.e., planar anchoring) without an energy

penalty.

The attractive part of the fluid-substrate potential is mod-

eled via a Yukawa potential. This has the advantage that by

204702-12 M. Greschek and M. Schoen J. Chem. Phys. 135, 204702 (2011)

varying the screening parameter, we are able to gradually

change the range of the fluid-substrate attraction. To compare

data for various screening parameters, we introduce two ad-

ditional parameters that allow us to fix both the position of

the minimum of the fluid-substrate potential and the value of

the potential at the position of that minimum. For purely tech-

nical reasons, we “confine” the liquid crystal to a slit-pore.

However, we make the distance between the pore walls suffi-

ciently large so that molecules do not interact simultaneously

with both walls. In effect, this permits us to study orienta-

tion effects in liquid-crystalline materials wetting single solid

surfaces.

We consider a single temperature T=0.9 which is both

below the bulk gas-(isotropic) liquid critical point and above

the gas-isotropic-nematic triple point. The latter is inferred

from the fact that with increasing μ, the transition to the

isotropic liquid and eventually from that phase to the nematic

are well separated in the bulk. This is concluded from plots in

Figs. 3and 4, which show that both transitions occur at about

the same μirrespective of the specific anchoring scenario. In

other words, the location of these transitions is unaffected by

the presence of the substrates.

Our key findings can be summarized as follows. Sub-

strates with degenerate anchoring are always covered with a

submolecularly thin film at sufficiently low μ. Increasing μ

further causes the thickness of this film to increase continu-

ously until a thick film consisting of a few molecular layers

has formed. If anchoring is degenerate but results in planar

orientation of the adsorbed molecules, a prewetting transition

occurs at which a submonolayer film is transformed discon-

tinuously into a molecularly thick film. A prewetting transi-

tion is also observed for the two monostable (d and h) anchor-

ing scenarios considered in this study. In this latter case, μij

is roughly the same for both d and h anchoring and exceeds

μij

xfor p anchoring significantly.

Anchoring is associated with orientational configura-

tional entropy. This can be realized by comparing degener-

ate with monostable anchoring. Specifically, for n anchoring,

where no orientation of a molecule with respect to the sub-

strate plane receives any energy penalty, more configurations

of molecules are possible than for both d and h anchoring.

Because no molecular orientation is preferred by n anchor-

ing, one might also consider this anchoring scenario as weak.

In these latter cases, all orientations except for a single one

receive such an energy penalty to a larger or lesser extend de-

pending on the degree of misalignment of molecules with the

easy axis. Moreover, comparing n with p anchoring, orienta-

tional configurational entropy should be larger for the former

because for the latter anchoring scenario only those orienta-

tions that lie on the unit circle in the x–yplane avoid receiv-

ing an energy penalty. Orientational configurational entropy is

therefore expected to decrease in the direction n →p→d, h.

Because no prewetting occurs at η=0.25 for n anchoring, a

prewetting transition is observed for p anchoring, and because

μij

xis smaller for p compared with both d and h anchoring, we

conclude that prewetting is shifted to lower μ, the larger the

orientational configurational entropy is. This interpretation is

further corroborated by the observation that μij

xis roughly the

same for the two monostable anchoring scenarios which are

expected to give rise to a similar orientational configurational

entropy. Hence, because of the close relationship between ori-

entational effects caused by different anchoring and the ex-

istence/location of the prewetting transition, we refer to this

phase transition as “orientational prewetting.” Hence, we con-

clude that orientational prewetting is caused by substrates re-

ducing the orientational configuration entropy of an adjacent

sufficiently thin liquid-crystalline film. As μincreases beyond

its value at the prewetting transition, the films at both sub-

strates (mono- and multistable anchoring) continue to grow

in thickness until finally complete wetting sets in.

Whether or not a thermodynamically stable prewetting

film forms depends on the range of the fluid-substrate attrac-

tion, i.e., whether one is above or below the prewetting crit-

ical temperature which is expected to depend on the fluid-

substrate attraction. These findings are qualitatively similar

to the ones reported earlier by Telo da Gama using density-

functional theory at mean-field level.20 Here we find that for

the shortest-range attraction (η=1.00), orientational prewet-

ting does not occur. Instead, the monolayer film adsorbed on

the monostable anchoring substrate vanishes during a wet-

ting transition when the entire liquid crystal between both

substrates condenses (complete wetting). The new liquid-like

phase is mostly isotropic except for the first few layers ad-

sorbed on the substrate which retain a certain degree of lo-

cal nematic order. This situation may be perceived as one in

which a substrate-controlled nematic region is completely wet

by an isotropic liquid. The decay of nematic order with in-

creasing distance from the substrate is closely correlated with

the range of fluid-solid attraction.

At a sufficiently high μ, the entire liquid crystal under-

goes an IN phase transition. The location of the IN phase tran-

sition is nearly independent of the precise anchoring scenario

as one would expect if the distance between the solid sub-

strates is sufficiently large as we have argued above. However,

both “snapshots” of characteristic configurations of molecules

as well as properly defined orientation profiles show that the

specific anchoring scenario has a significant impact on 

Hence, even though the range of direct fluid-substrate inter-

action is finite, the specific anchoring scenario is imprinted

on more remote portions of the liquid crystal (i.e., portions

farther removed from the solid substrate) through mediation

by molecules that are closer to the substrate.

ACKNOWLEDGMENTS

We are grateful for financial support from the Inter-

national Graduate Research Training Group 1524 “Self-

Assembled Soft Matter Nanostructures at Interfaces.”

APPENDIX: PROPERTIES OF THE

GRAND-POTENTIAL DENSITY

Here we introduce additional properties of the grand-

potential density that are useful in analyzing the phase behav-

ior of fluid phases. First, one notices from Eqs. (3.4) and (3.6)

that ∂ωi

∂μ {·}\μ

=−ρi<0,(A1)

204702-13 Orientational prewetting by a liquid crystal J. Chem. Phys. 135, 204702 (2011)

where “{ ·}\μ” is the shorthand notation to indicate that upon

differentiation all other natural variables of ωexcept for μ

are to be held constant. Subscript irefers to the specific ther-

modynamic phase under consideration. Because the number

density is a positive definite quantity by definition, we an-

ticipate plots of ωiversus μat constant T,A, and szto be

monotonically decreasing where the slope is the larger the

higher ρiis for a given phase. Second, using the definition

κ≡−A−1(∂A/∂P)T,N,s

zfor the transverse isothermal com-

pressibility together with the Gibbs-Duhem equation40

0=SdT+Ndμ−Asz0dP+PzzA0dsz,(A2)

one can easily verify that at constant Tand sz40

∂2ωi

∂μ2{·}\μ

=−ρ2

iκ<0,(A3)

which holds because thermodynamic stability requires κ

>0.58 Hence, for T=constant, the functions ωi(μ) are con-

cave, that is the inequality

ωi(μ)≥ωi(μ)(μ−μ)+ωi(μ)(μ −μ)

μ −μ,(A4)

holds for all μ≤μ≤μ. The choice of the interval [μ,μ]

is arbitrary but must belong to the range of chemical poten-

tials for which phase iis at least metastable. The important

conclusion one may draw from Eqs. (A1) and (A3) is that if

the equation

ωiμij

x=ωjμij

x,T=constant, (A5)

has a solution, there will be one and only one solution μij

which is the chemical potential at which phases iand jcoex-

ist. Moreover, for subcritical temperatures, we anticipate ρix

= ρjxin general. If ρix<ρ

jx(ρix>ρ

jx), phase iwill be the

more (less) stable one for μ<μ

xwhereas phase jwill be the

more (less) stable one for μ>μ

1B. Jerôme, Rep. Prog. Phys. 54, 391 (1991).

2P. Sheng, Phys. Rev. Lett. 37, 1059 (1976).

3P. Sheng, Phys. Rev. A 26, 1610 (1982).

4K. Miyano, Phys. Rev. Lett. 43, 51 (1979).

5E. Grelet and H. Bock, Europhys. Lett. 73, 712 (2006).

6A. A. Sonin, The Surface Physics of Liquid Crystals (Gordon and Breach,

Amsterdam, 1995).

7Y. Yi, G. Lombardo, N. Ashby, R. Barberi, J. E. Maclennan, and

N. A. Clark, Phys. Rev. E 79, 041701 (2009).

8H. Hsiung, T. Rasing, and Y. R. Shen, Phys. Rev. Lett. 57, 3065 (1986).

9M. P. Valignat, N. Fraysse, and A. M. Cazabat, Langmuir 9, 3255 (1993).

10W. Chen, L. J. Martinez-Miranda, H. Hsiung, and Y. R. Shen, Phys. Rev.

Lett. 62, 1860 (1989).

11R. Lucht and C. Bahr, Phys. Rev. Lett. 80, 3783 (1998).

12C. Bahr, Europhys. Lett. 88, 46001 (2009).

13X. Liu, D. W. Allender, and D. Finotello, Europhys. Lett. 59, 848 (2002).

14A. Poniewierski and T. J. Sluckin, Mol. Cryst. Liq. Cryst. 111, 373 (1983).

15T. J. Sluckin and A. Poniewierski, Phys. Rev. Lett. 55, 2907 (1985).

16A. Mauger, G. Zribi, D. L. Mills, and J. Toner, Phys. Rev. Lett. 53, 2485

(1984).

17A. K. Sen and D. E. Sullivan, Phys.Rev.A35, 1391 (1987).

18H. Stark, J.-I. Fukuda, and H. Yokoyama, J. Phys.: Condens. Matter 16,

S1911 (2004).

19P. Patrício, C.-T. Pham, and J. J. Romero-Enrique, Eur. Phys. J. E 26,97

(2008).

20M. M. Telo da Gama, Mol. Phys. 52, 611 (1984).

21L. Onsager, Ann. N.Y. Acad. Sci. 51, 627 (1949).

22D. de las Heras, L. Mederos, and E. Velasco, Phys.Rev.E68, 031709

(2003).

23D. de las Heras, E. Velasco, and L. Mederos, J. Chem. Phys. 120, 4949

(2004).

24M. Dijkstra, R. van Roij, and R. Evans, Phys. Rev. E 63, 051703 (2001).

25R. van Roij, M. Dijkstra, and R. Evans, Europhys. Lett. 49, 350 (2000).

26R. Zwanzig, J. Chem. Phys. 39, 1714 (1979).

27M. P. Allen, J. Chem. Phys. 112, 5447 (2000).

28M. T. Downton and M. P. Allen, Europhys. Lett. 65, 48 (2004).

29L. Bellier-Castella, D. Caprion, and J.-P. Ryckaert, J. Chem. Phys. 121,

4874 (2004).

30E. A. Müller, J. Rodriguez-Ponce, A. Qualid, J. M. Romero-Enrique, and

L. F. Rull, Mol. Simul. 29, 385 (2003).

31F. Barmes and D. J. Cleaver, Phys. Rev. E 69, 061705 (2004).

32D. L. Cheung, J. Chem. Phys. 128, 194902 (2008).

33Z. Zhang, A. Chakrabarti, O. J. Mouritsen, and M. J. Zuckermann, Phys.

Rev. E 53, 2461 (1996).

34S. Hess and B. Su, Z. Naturforsch. 54a, 559 (1999).

35M. Greschek, M. Melle, and M. Schoen, Soft Matter 6, 1898 (2010).

36M. Greschek and M. Schoen, Soft Matter 6, 4931 (2010).

37M. Greschek and M. Schoen, Phys.Rev.E83, 011704 (2011).

38H. Steuer, S. Hess, and M. Schoen, Physica A 328, 322 (2003).

39H. Steuer, S. Hess, and M. Schoen, Phys.Rev.E69, 031708 (2004).

40M. Schoen and S. H. L. Klapp, Reviews in Computational Chemistry

(Wiley-VCH, Hoboken, 2007), Vol. 24.

41M. Schoen, O. Paris, G. Günther, D. Müter, J. Prass, and P. Fratzl, Phys.

Chem. Chem. Phys. 12, 11267 (2010).

42C. G. Gray and K. E. Gubbins, Theory of Molecular Fluids (Clarendon,

Oxford, 1984), Vol. 1, Chap. E 2.

43I. Pardowitz and S. Hess, Physica A 100, 540 (1980).

44W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numeri-

cal Recipes in FORTRAN (Cambridge University Press, Cambridge, 1989),

Chap. 11.1.

45A. Richter and T. Gruhn, J. Chem. Phys. 125, 064908 (2006).

46P. G. de Gennes and J. Prost, The Physics of Liquid Crystals (Oxford Sci-

ence Publications, Oxford, 1995).

47R. Eppenga and D. Frenkel, Mol. Phys. 52, 1303 (1984).

48W. Maier and A. Saupe, Z. Naturforsch 14a, 882 (1959).

49W. Maier and A. Saupe, Z. Naturforsch. 15a, 287 (1960).

50G. Arfken, Mathematical Methods for Physicists (Academic, San Diego,

1985).

51R. Low, Eur. J. Phys. 23, 111 (2002).

52T. Gruhn and M. Schoen, Phys.Rev.E55, 2861 (1997).

53D. J. Adams, Mol. Phys. 29, 307 (1975).

54M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids

(Clarendon, Oxford, 1987).

55J. H. Thurtell, M. M. Telo da Gama, and K. E. Gubbins, Mol. Phys. 54, 321

(1985).

56E. M. del Rio, M. M. Telo da Gama, and E. de Miguel, Phys. Rev. E 52,

5028 (1995).

57L. F. Rull, J. M. Romero-Enrique, and E. A. Müller, J. Phys. Chem. C 111,

15998 (2007).

58H. B. Callen, Thermodynamics (Wiley, New York, 1960).