A Equilibrium grand–canonical simulations for the monolayer in the lattice model

Monolayers of hard rods on planar substrates: I. Equilibrium

Abstract

The equilibrium properties of hard rod monolayers are investigated in a lattice model (where position and orientation of a rod are restricted to discrete values) as well as in an off–lattice model featuring spherocylinders with continuous positional and orientational degrees of freedom. Both models are treated using density functional theory and Monte Carlo simulations. Upon increasing the density of rods in the monolayer, there is a continuous ordering of the rods along the monolayer normal (“standing up” transition). The continuous transition also persists in the case of an external potential which favors flat–lying rods in the monolayer. This behavior is found in both the lattice and the continuum model. For the lattice model, we find very good agreement between the results from the specific DFT used (lattice fundamental measure theory) and simulations. The properties of lattice fundamental measure theory are further illustrated by the phase diagrams of bulk hard rods in two and three dimensions.

pacs:

I Introduction

Several systems of scientific and technological interest can be characterized as being monolayers of anisotropic particles, such as Langmuir monolayers (1) or very thin films of elongated organic molecules (2) such as organic semiconductors (3); (4). Since the particular molecular interactions in these systems may be very complicated, it is worthwhile to investigate simpler models of anisotropic colloids to obtain general insights into the thermal behavior of these systems. Among these, hard–body models (where particles interact only via their excluded volume) are a natural starting point to assess effects of anisotropy, both for thermal equilibrium and non–equilibrium conditions (i.e. growth of the monolayer).

We present our investigations on equilibrium and on growth of hard–rod monolayers in two papers, where in the present first one, we focus on equilibrium properties and in the second, we treat the growth process. In both papers, the focus will be on hard–rod lattice models since for these a fairly transparent theoretical analysis of equilibrium and growth is possible in the framework of density functional theory (DFT). In particular, we use the framework provided by fundamental measure theory (FMT) (5) within which very accurate density functionals for systems of anisotropic hard particles have been constructed for continuous (6) and, important for the present investigations, also for lattice models (7); (8). For the subsequent investigations of the growth process, kinetic Monte–Carlo simulations on a lattice are a natural starting point. Since lattice models inevitably restrict the translational and orientational degrees of freedom of rods, we will also present results for an off–lattice hard rod (spherocylinder) model and identify similarities and differences between lattice and off–lattice models.

The restriction of orientation in hard–rod (cuboid) models comes with the benefit that density functionals from FMT become tractable and therefore also analytic results can be derived. For continuous translational degrees of freedom and in three dimensions (3d) a rich phase diagram was derived (9). Although not all details are the same in a simulated phase diagram of hard cuboids with unrestricted orientation (10), the restricted orientation model gives a good first estimate of what can be expected. If the particles are restricted to a plane (the monolayer case), the first–order isotropic–nematic transition becomes a continuous one, according to FMT in the restricted orientation model (11); (12). The orientational order perpendicular to the plane is proportional to the density for low densities. An approximate DFT and simulations for hard ellipsoids (unrestricted orientation) seem to confirm this behavior although very low densities have not been sampled in the simulations (13). Such a possible qualitative change of the nature of the nematic transition through dimensional restriction is very interesting by itself, and therefore we will establish analytically the behavior explicitly in the low–density limit for both lattice and continuum models. The presence of an orientation–dependent external potential in the monolayer plane (substrate potential) does not change the continuous nature of the nematic transition but the onset of particles “standing up” may become very sharp for substrate potentials which actually favor particles “lying down”.

The structure of the paper is as follows: In Sec. II, we describe the lattice version of fundamental measure theory (FMT) for hard rod mixtures and give illustrative examples for the functionals. The bulk equilibrium properties of monocomponent rods in two dimensions (2D) and three dimensions (3D) are briefly discussed, followed by the results for the monolayer (3D confined). Sec. III discusses the spherocylinder off–lattice model for the monolayer using DFT in the low–density limit and using simulations. Sec. IV discusses similarities and differences between the lattice and off–lattice models and gives a summary. Two appendices briefly discuss the grand canonical simulation method for the lattice model and the derivation of the excluded area between hard rods (in the continuum model) whose centers are confined to a plane.

Ii Density Functional Theory for hard rod lattice models

ii.1 Fundamental measure theory

The rod model used in this work is formulated on a simple cubic lattice in dimensions. A lattice point is specified by a set of integers (). The lattice constant is the unit of length. Hard rods are lines (1D), rectangles (2D) or parallelepipeds (3D) with corners sitting on lattice points and thus their geometry is specified by their extent in the cartesian directions which are again sets of integers. The position of a rod is specified by the corner whose lattice coordinates are minimal each (see Fig. 1). Hard rods are not allowed to overlap (but they may “touch”, i.e. share surfaces), thus the interaction potential for two rods and of species and at positions and with extensions and is given by

(1)

Here, is the rod overlap function given by

(2)
(3)

The overlap function is 1 whenever there is overlap in all lattice dimensions, meaning that the rods are disjunct for (see. Fig. 1). Note that due to the chosen convention for the rod location the overlap function is not symmetric in the rod locations and .

Figure 1: (color online) Definitions for the example of hard 2x3–rods in . Rod location is specified by the position of the lower left corner (i.e., the corner whose lattice coordinates are minimal). Rods may “touch” (left picture) but not overlap (right).

In the following, we consider such a rod mixture with species subject to external fields where acts on rod species . At lattice site , the number density of rods per lattice site is specified by where is the density of rod species , i.e. the probability of a given site to be occupied by the lower left corner of a particle. In density functional theory, all equilibrium properties of a rod mixture in external fields are obtained by minimizing the grand potential functional

(4)

with respect to the particle densities . The chemical potential for rod species is denoted by If different species belong to the same type of rod in different orientations, the corresponding chemical potentials must be equal in equilibrium. denotes the ideal gas contribution to the free energy functional, given by

(5)

Energies are measured in units of throughout the paper.

The exact form of the excess free energy functional is in general unknown, in this work we will approximate it within the fundamental measure approach. For lattice models of hard rods, this approach has been worked out in Refs. (7); (8), resulting in an approximative form for which we apply in the present study (Lafuente–Cuesta functional).

Figure 2: The four FMT weight functions for a rod with edge lengths . The lattice point at which the weight functions are evaluated is denoted by . The thick points indicate on which lattice points the weight function is 1.

The class of free energy functionals derived in Refs. (7); (8) makes use of weighted densities which are defined as convolutions of densities with weight functions :

(6)

Convolutions () on the lattice are defined as

(7)

The –dimensional index specifies different weight functions , with allowed values only. The weight functions (specific for species ) have the meaning of defining a support of rods with edge lengths , i.e. they are 1 on points covered by and 0 otherwise. This can be formalized using the –function already employed for defining rod overlap (Eq. (3)),

(8)

The edge lengths of rods are related to those of the rods as follows:

(9)

i.e. whenever the index is 0, the edge length of in the dimension is shortened by 1 compared to the corresponding edge length of , otherwise () the edge length is identical. In particular, for all rods are identical to . The corresponding weighted densitiy has the meaning of a local packing or volume fraction of rods at point . Fig. 2 illustrates the four possible weight functions for a rods with edge lengths on a 2D lattice.

As a second ingredient, the Lafuente–Cuesta functional needs the excess free energy of a zero–dimensional (0d) cavity, , i.e. a restricted domain on the lattice which can only hold one particle at a time. Such a cavity may consist of more than one point where the rod is positioned. Furthermore, for a mixture the set of points specifying the allowed location of species does not need to coincide with the corresponding set for species . Note that the sets and are not independent since the 0D cavity property is required to hold globally for the mixture and not just for the individual components. The free energy of such a cavity is a function only of the total packing fraction in the cavity,

(10)

Using this 0d free energy, the Lafuente–Cuesta excess free energy functional is given by

(11)

Remember that is a –dimensional index with entries only. In Eq. (11), and is the difference operator whose action on a function is given by .

It can be shown that as defined above yields the correct excess free energy, Eq. (10), for any 0D cavity (7); (8). In order to assess the accuracy of the expression for situations of less severe confinement, we evaluate explicitly the properties of different bulk systems in Sec. II.3. In a first step, however, we illustrate the construction of the Lafuente–Cuesta functional by applying it to different mixtures in 1D, 2D, and 3D.

ii.2 Special cases

Here we give the explicit functionals for some special mixtures. The equilibrium properties of examples (b) and (c) (2D and 3D systems) will be discussed in Sec. II.3 and those of example (d) (monolayer) in Sec. II.4.

  • : Mixture of hard rods in one dimension. The excess free energy functional is given by

    (12)

    This is the well–known exact solution for the 1d lattice hard rod mixture, derived in Ref. (7) following the recipe from Ref. (14) which treats the 1D continuum hard rod mixture. Another yet different derivation can be found in Ref. (15).

  • : A system of rods with length and width 1 corresponds to the binary mixture with rod lengths and . The excess free energy functional is given by

    (13)

    The weighted densities are given by

    (14)

    Note that the weights since they correspond to the support of rods with width 0. Likewise .

  • : A system of rods with length and height/width 1 corresponds to the ternary mixture with rod lengths , and . The excess free energy functional is given by

    (15)

    The weighted densities are given by

    (16)
    (17)

    Similarly to case (b), the weights whenever and since they correspond to the support of rods with width 0.

  • (confined), the monolayer: A system of rods with length and height/width 1 whose positions are constrained to a 2D–plane corresponds to a 2D ternary mixture with rod lengths , (rods lying in–plane) and (rods standing up). The excess free energy functional is given by formally the same functional as in (a),

    (18)

    but now the weighted densities are given by

    (19)

ii.3 Equilibrium bulk properties in 2D and 3D

Figure 3: (color online) (a) Rods in : Demixing order parameter as a function of the total packing fraction for different rod lengths . Dotted lines correspond to the approximate solution near the onset of demxing (Eq. (26)).

, the binary mixture with rod lengths and

In the bulk, both densities ( and ) and all weighted densities are constant. We introduce the total density and denote by the total packing fraction. Furthermore , and is an order parameter for the demixed state. We refrain from calling a nematic order parameter, since the alignment of rods corresponds just to a demixed state between species 1 and 2, and the corresponding transition has the character of a liquid–vapor transition (16). The bulk free energy density, , written in dependence on the variables and becomes

(20)
(21)
(22)
(23)

At fixed , the equilibrium demixing parameter is found by solving . For , the mixed state () is the only solution and is minimal there, For there exists a critical packing fraction above which three solutions signal demixing: the solutions have lower free energy. At , there is no jump in the demixing parameter which is the behavior also observed at a liquid–vapor transition. Therefore one may expand

(24)

and find the critical packing fraction by solving , with the solution

(25)

The equilibrium demixing in the vicinity of can be approximated by solving for using the Taylor approximation (24), giving

(26)

The behavior of near born out by the approximate theory is of course of mean-field type.

These findings can be compared with simulation work which finds the demixing transition for (17) and a critical packing fraction (19). Thus, FMT overestimates the tendency to demix. Note, however, that the demixing follows from a single functional, unlike other approaches which assume distinct epxressions for the isotropic and the demixed phase free energies (19).

, the ternary mixture with rod lengths , and

The total density is and the total packing fraction is . We define the order parameters

(27)

signifies an excess () or depletion () of particles in –direction (nematic state) while signals order in the –plane orthogonal to the nematic director (biaxial state). The bulk free energy density, , written in dependence on the variables , and becomes

(28)
(29)
(30)
(31)
(32)

Minimization of the total free energy density with respect to and shows that the model has a stable nematic state (, ) for . Note that the director could also be oriented along the – or –axis instead of the chosen –axis. A pure nematic state with director along the –axis and order parameter is equivalent to a minimum free energy state with and using the order parameters (27). Therefore this is not a biaxial state. The associated liquid–nematic transition is of first order, and we have determined coexistence between the liquid and the nematic state by performing the common tangent construction for the free energy density (liquid phase) and (nematic phase) which implies equality of the chemical potential and pressure . Results are shown in Fig. 3(b). The packing fractions of the coexisting nematic state are very well described by , and the gap in packing fractions of the coexisting states has a maximum of at and tends to zero as .

Figure 4: (color online) (a) Order parameter for rods standing up vs. total density. Lines are results from FMT and symbols are results from Monte Carlo simulations reported in Ref. (20). Thin lines are results from our GCMC simulations where a running average of 20 points on density intervals of 0.04 has been taken. (b) Phase diagram from FMT showing a reentrant behavior for mixing () and demixing () in the plane. The rod length is treated as a continuous variable. The critical point occurs for a rod length of at a density of .

ii.4 The monolayer system

In the lattice model, this is the effectively 2D ternary mixture with rod lengths , (rods lying in–plane) and (rods standing up).

The total density is and the total packing fraction in the plane is . The order parameters and are the same as in Eqs. (27). signifies an excess of particles “standing–up” (nematic state) while signals demixing of “lying–down” particles (biaxial state, if additionally ). In the bulk free energy density, , one can identify whereas the excess part becomes

(33)

At fixed total density , the minimization of the free energy with respect to and gives the following picture. For “small” rod lengths there is no biaxial state (, no demixing in the plane) but the “nematic” order parameter monotonically and smoothly grows from 0 to 1 when the total density varies between 0 and 1 (close-packed state of rods standing up). Results for are shown in Fig. 4(a) which demonstrates that for increasing the rods quickly “stand up”. The FMT results show excellent agreement with Monte Carlo simulation results (20) on the same confined model for and 6. For larger rod lengths ( and 10) the agreement with our grand canonical Monte Carlo (GCMC) simulations is only slightly worse. The implementation of GCMC is briefly described in App. A.

For , FMT predicts reentrant demixing in the plane, i.e. in a certain interval for the total density the biaxiality parameter will be nonzero, . This reentrant behavior is qualitatively understood as follows. In the model it was found that the critical density of demixing of planar rods is . For increasing one therefore expects . On the other hand, for a certain but increasing the fraction of planar rods initially grows, reaches a maximum and becomes smaller again since the rods stand up, see Fig. 4(a). Therefore, if there exists a lower demixing density then one would expect the existence of a higher remixing density owing to the reduction of the planar rod density. As in the model, the demixing transition is continuous and therefore the densities can be found by the following argument: Let and be chemical potentials for the order parameters and . For a mixed state (), we define through . As before, we may expand

(34)

At the de-/remixing densities one has the condition

(35)

which needs to be solved numerically. The results are shown in Fig. 4(b), showing the onset of demixing at and a maximum density interval for the demixed state at around .

The continuous behavior of and the reentrant demixing are actually very similar to the behavior found in the FMT study of the restricted orientation model with continuous translational degrees of freedom (11). There biaxial ordering sets in at larger rod lengths, .

The results in Fig. 4 suggest , i.e. the continuous nematic ordering sets in at . This is easily understood in a low–density expansion of the FMT excess free energy (33) which is exact up to second order. Assuming no biaxiality () and combining ideal and excess part we find for the free energy derivative with respect to :

(36)

Note that in the excess part of , at fixed density, there is a constant term driving the system to for . This is different from the 2D and 3D bulk systems where this constant term is absent and thus (for low densities) is always unfavorable in terms of free energy cost. The equilibrium solution at is found as

(37)

Hence, for large the lattice model predicts a scaling .

Figure 5: (color online) Order parameter for rods standing up vs. total density subject to a substrate potential (rod length ). The substrate potential is parametrized as per unit length such that , and thus . Lines are DFT results, symbols results from GCMC simulations (see App. A). The error is smaller than the symbol size.

Finite substrate potential

One may ask whether a finite substrate potential could alter the continuous transition found above. It is natural to assume that the substrate potential acts equally on the flat–lying species 1 and 2 and differently on the upright species 3. Hence the external contribution to the free energy becomes

(38)

Therefore the free energy derivative with respect to is modified as with . For the ideal gas limit this implies an initial ordering on the substrate with order parameter

(39)

If the substrate is strongly attractive for the flat–lying species 1 and 2 (), then we find . At nonzero densities, the solution of (Eq. 36 with the external contribution) is obtained in the form . For small deviations from equilibrium, , we can invert this function and obtain

(40)

with and . Although the range of validity is very limited, it implies that the qualitative behavior for is unchanged since the slope of is always positive. Thus the transition stays continuous. However, for increasing the “standing up” transition of the monolayer becomes increasingly steep at moderate densities, see Fig. 5 where we show the behavior for . For these moderate densities the expansion up to second order is not valid anymore. Especially for the case the behavior near looks as if has a bifurcation point here, similar to the demixing transition in the 2D bulk system discussed in Sec. II.3. The density at which this apparent transition occurs is the close–packing density for rods lying flat. However, for finite potentials it is not a phase transition since maintains its linear behavior of with nonzero slope at very small densities.

Iii Monolayers of hard spherocylinders

For the lattice monolayer discussed in the previous section it does not matter which rod point or segment is actually fixed to the plane since all choices lead to the same effective 2D model. Physically, fixing the end point corresponds to the case of rods on a hard substrate while fixing some point in the middle of the rod applies to Langmuir monolayers. For a continuous model of hard rods, there should be a difference between the two cases which is not expected to be qualitative (with regard to the type of transition). As can be seen below, the low–density behavior of long rods with large aspect ratios is actually insensitive to the choice of confining plane. Therefore we will present simulation results below only for the case of fixed mid–points.

iii.1 DFT in an expansion up to second order in density

We consider hard sperocylinders with length and diameter whose centers or ends are fixed on a plane. In order to investigate the nature of the orientation transition, we consider a low–density expansion of the free energy. For the well–studied model of hard rods in 3D, this method was used to establish the onset of nematic order as a bifurcation and the nature as a first order transition (21). The free energy density up to second order in density, including the contribution from an external potential, is given by

(41)
(42)
(43)
(44)

Here, is an inhomogeneous particle density in two dimensions and units of [length] which depends on the space point and the orientation of the rod , specified by the polar angle and the azimuthal angle . The integral over orientations is defined as

(45)

is the thermal de–Broglie length. is the overlap function between rods for given orientations of and distance between the particles. It is 1 if there is overlap, otherwise zero. The external (substrate) potential is measured in units of .

We consider only orientation–dependent substrate potentials, , and bulk states, i.e. no space dependence of the density and introduce the orientation distribution :

(46)

Then the ideal, excess and external part of the free energy per particle () become:

(47)
(48)
(49)

Here, is the excluded area between the rod centers (or ends) with fixed orientations of the rods.

In equilibrium, minimizes . From we obtain

(50)

where is a constant which must ensure that is properly normalized, i.e. . It is determined by exponentiating Eq. (50) and integrating over :

(51)

The orientation–dependent substrate potential gives rise to a non–constant orientation distribution in the ideal gas limit:

(52)

which is normalized to 1. We introduce the small deviation and linearize Eq. (51) in :

(53)

The constant ensures the necessary normalization condition . If one expands in powers of then one finds the leading order solution

(54)

This expression is equivalent to Eq. (40) in the lattice model and shows that any deviations from the ideal gas distribution are continuous and proportional to the density .

In the absence of a substrate potential (), we can proceed further. Without loss of generality, we put rod 1 at the coordinate center with orientation (director) . Rod 2 has the director . The excluded area depends in general on the three angles . If we consider only nematic order without biaxiality, , then we can define the integrated overlap area

(55)

If we take the polar angle (with respect to the interface normal) in the interval , symmetry considerations give us . Since also , the integration domain over can be restricted to . The nematic order parameter in the monolayer is defined by

(56)

where is the second of the Legendre polynomials . It is also useful to introduce the Legendre coefficients of the excluded area:

(57)

Owing to the symmetry of the excluded area, only projections onto even Legendre polynomials are nonzero. Using these definitions, the nematic order parameter in the case of no substrate potential is obtained by projecting with onto the solution for in Eq. (54):

(58)

This is an interesting result since it tells us that as long as the leading off–diagonal Legendre coefficient of the excluded area is nonzero. This is precisely the case in the monolayer system (see below), whereas in 3D this coefficient vanishes. The linearity is completely equivalent to the linearity found in the lattice model in the absence of a substrate potential (see Eq. (37)).

The linearized equation (53) is connected to an approximated free energy per particle through . is quadratic in and is defined to give the the difference to the isotropic state:

(59)

For the leading order solution (54), the free energy can be explicitly evaluated. It is convenient to use the Legendre expansion of the solution: with . Using furthermore