Generalized Onsager theory for strongly anisometric patchy colloids H. H. Wensink and E. Trizac Citation: The Journal of Chemical Physics 140, 024901 (2014); doi: 10.1063/1.4851217 View online: http://dx.doi.org/10.1063/1.4851217 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/2?ver=pdfcov Published by the AIP Publishing

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

THE JOURNAL OF CHEMICAL PHYSICS 140, 024901 (2014)

Generalized Onsager theory for strongly anisometric patchy colloids H. H. Wensink1,a) and E. Trizac2 1

Laboratoire de Physique des Solides - UMR 8502, Université Paris-Sud and CNRS, 91405 Orsay Cedex, France 2 Laboratoire de Physique Théorique et Modèles Statistiques - UMR 8626, Université Paris-Sud and CNRS, 91405 Orsay Cedex, France

(Received 3 September 2013; accepted 5 December 2013; published online 8 January 2014) The implications of soft “patchy” interactions on the orientational disorder-order transition of strongly elongated colloidal rods and flat disks is studied within a simple Onsager-van der Waals density functional theory. The theory provides a generic framework for studying the liquid crystal phase behaviour of highly anisometric cylindrical colloids which carry a distinct geometrical pattern of repulsive or attractive soft interactions localized on the particle surface. In this paper, we apply our theory to the case of charged rods and disks for which the local electrostatic interactions can be described by a screened-Coulomb potential. We consider infinitely thin rod like cylinders with a uniform line charge and infinitely thin discotic cylinders with several distinctly different surface charge patterns. Irrespective of the backbone shape the isotropic-nematic phase diagrams of charged colloids feature a generic destabilization of nematic order at low ionic strength, a dramatic narrowing of the biphasic density region, and a reentrant phenomenon upon reducing the electrostatic screening. The low screening regime is characterized by a complete suppression of nematic order in favor of positionally ordered liquid crystal phases. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4851217] I. INTRODUCTION

Many colloidal dispersions, such as natural clays, and (bio-)macromolecular systems consist of rod- or disk-shaped mesogens whose intrinsic ability to form liquid crystalline order gives rise to unique rheological and optical properties.1 Despite their abundance in nature, the statistical mechanics of fluids containing non-anisometric particles in general (and oblate ones in particular) has received far less attention than that of their spherical counterparts. The possibility of a firstorder disorder-order transition from an isotropic to a nematic phase was first established theoretically by Onsager2 in the late 1940s. Although originally devised for rod-like particles in solution, his theory also makes qualitative predictions for plate-like particles based on the central idea that orientationdependent harshly repulsive interactions alone are responsible for stabilizing nematic order. Subsequent numerical studies have fully established the phase diagram of hard prolate3–8 and oblate hard cylinders.9–11 Owing to the simplicity of the interaction potential hard-body systems constitute an essential benchmark for the study of liquid crystals and their phase stability. Temperature becomes merely an irrelevant scaling factor in the free energy and the phase behaviour is fully determined by the volume fraction occupied by the particles and the aspect ratio. At high volume fraction, additional entropydriven disorder-order transitions occur where a nematic fluid transforms into positionally ordered phases.12 Depending on the cylinder aspect ratio the system may develop a smectic phase, characterized by a one-dimensional periodic modulation along the nematic director, or a columnar phase consisting of columns with a liquid internal structure self-assembled a) [email protected]

0021-9606/2014/140(2)/024901/13/$30.00

into a two-dimensional crystal lattice. Similar to nematic order, the formation of smectic, columnar, or fully crystalline structures is based entirely on entropic grounds;13 the loss of configurational entropy associated with (partial) crystalline arrangement is more than offset by a simultaneous increase in translational entropy, that is, the average free space each particle can explore becomes larger in the ordered phase. In most practical cases, however, particle interactions are never truly hard and additional enthalpic contributions play a role in the free energy of the system. Long-ranged interactions usually originate from the presence of surface charges leading to electrostatic repulsions between colloids14, 15 or traces of other colloidal components such as non-adsorbing polymers, which act as depletion agents and give rise to effective attractive interactions.16, 17 Other site-specific interactions may originate from hydrogen-bonding18 or endfunctionalized polymers such as DNA grafted onto the colloid surface.19 Depending on their nature (repulsive or attractive), interaction range, and topological arrangement on the particle surface, these site-specific directional interactions may greatly affect the self-assembly properties of anisometric particles.20–22 In this context, it is also worth mentioning recent progress in the fabrication of anisometric colloids with “patchy” interactions23, 24 where the interplay between patchiness and the anisometric backbone shape offers a rich and intriguing repertoire of novel structures.25 These recent developments suggest the need for a comprehensive theory for lyotropic systems which explicitly accounts for these patchy interactions. The aim of the present paper is to set up such a theory by combining the classic Onsager theory for slender hard bodies with a meanfield van der Waals treatment for the additional longranged interactions.26–30 Most molecular-field type theories

140, 024901-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

024901-2

H. H. Wensink and E. Trizac

developed to date focus on rod-like mesogens with dispersion interactions represented by an orientation-dependent potential with some radially symmetric spatial variation, akin to a Maier-Saupe form.31–34 Here, we shall lay out the framework for the more general case of slender rod and diskshaped cylinders carrying site interactions with arbitrary integrable form and spatial arrangement. By exploiting the simple second-virial structure of the Onsager reference free energy, we show that these soft patchy interactions, on the mean-field level, give rise to a non-trivial orientationdependent van der Waals (or molecular field) term which strongly affects the disorder-order transition in the fluid state. We illustrate its practical use by focusing on isotropicto-nematic and nematic-to-smectic or columnar phase transitions in systems of charged prolate and discotic colloid in the salt-dominated regime, a subject of considerable research interest given that natural clays consist of strongly charged colloids. The majority of clays are composed of sheet-like minerals colloids35, 36 but rod-shaped mineral colloids may display similar properties.37–39 It is still largely unclear how the interplay between particle shape and electrostatics controls the structure and dynamics of clay systems. The fundamental understanding is further complicated by the fact that both the magnitude and sign of the local charge density may vary significantly along the particle surface. For instance, under certain chemical conditions laponite platelets40 adopt opposite face and rim-charges and the intrinsic patchiness of the electrostatic interactions may lead to unusual liquid behaviour.41 Incorporating these patchy interactions into a state-of-the-art statistical physical machinery to extract structural information remains a daunting task. Headway can be made by using computer simulation in which context a number of coarse-grained models for non-isometric charged colloids have been studied over the past decade.42–47 With the present theory, we aim to set a first step towards linking microscopic patchiness of soft interactions to liquid crystal stability for strongly anisometric colloids. We apply the generalized Onsager theory to the case of charged cylinders interacting through an effective Yukawa potential and demonstrate a generic destabilization and non-monotonic narrowing of the biphasic gap upon reducing the electrostatic screening. The influence of the geometric pattern of the charge patches can be incorporated explicitly by means of a form factor as shown for disklike colloids. The present calculations, however, merely serve an illustrative purpose and the main goal is to open up viable routes to studying more complicated surface charge architectures of clay nano sheets48–51 or anisotropic Janus particles.52, 53 Moreover, the theory can be further refined by using effective parameters, pertaining to the backbone shape, charge density, screening constant, etc., in order to enable more quantitative predictions for highly charged anisometric colloids. Although the Onsager treatment is strictly limited to low to moderate density, it offers possibilities to assess the stability of high-density liquid crystal phases on the level of a simple bifurcation analysis.54 We show that it is possible to extend the generalized-Onsager form into a full density functional form by using a judiciously chosen parametric form for the one-body density. This holds promise for incorporat-

J. Chem. Phys. 140, 024901 (2014)

ing soft interactions into more sophisticated hard-body density functionals such as those based on fundamental measure theory,55–57 weighted-density approximations,8 renormalized Onsager theories,30, 58 or cell-theories.3, 7, 59 The use of reliable non-local reference free energy functionals is expected to give a more quantitative account of patchy rods or disks with broken translational symmetry induced by a high particle density, geometric confinement,60 or surfaces.61 The generalized Onsager theory bears some resemblance to other interaction-site models such as PRISM/RISM theories62, 63 which have been invoked to study the thermodynamic properties of isotropic plate fluids but have not yet proven capable of treating liquid crystal phases at higher particle densities. The effect of attractive interparticle forces on the bulk phase behaviour of ionic liquid crystals has been scrutinized in Ref. 64 using a mean-field theory of the Gay-Berne potential for ellipsoidal mesogens. The remainder of this paper is structured as follows. In Sec. II, we outline the mean-field Onsager theory for soft patchy cylinders with vanishing thickness. The theory will then be applied in Sec. III to study the isotropic-nematic phase diagram of charged rod- and disklike cylinders in the strong screening regime. Results for the isotropic-nematic phase diagrams will be presented in Sec. IV. Possible ways to include spatially inhomogeneous liquid crystals into the generalized Onsager treatment are highlighted in Sec. V. Finally, some concluding remarks are formulated in Sec. VI.

II. MEAN-FIELD ONSAGER THEORY FOR SOFT PATCHY POTENTIALS

Let us consider a system of N infinitely thin colloidal cylindrical disks or rods with length L and diameter D at positions {rN } and orientations {N } in a 3D volume V at temperature T. We assume the particle shape to be maximally anisotropic so that the aspect ratio L/D → ∞ (infinitely elongated rods) and L/D↓0 (infinitely flat disks). In the fluid state, the particle density ρ = N/V is homogeneous throughout space. Following Onsager’s classical theory,2 we may write the Helmholtz free energy as follows: βF ρ ∼ ln Vρ + ln 4πf () − N 2



 dr(r; 1 , 2 )

,

V

(1) with β −1 = kB T in terms of Boltzmann’s constant kB and V the total thermal volume of a cylinder including contributions from the rotational momenta. The brackets  denote orientation averages  ·  = df()( · ) and ·  = d1 d2 f (1 )f (2 )(·) in terms of the orientational distribution function (ODF) f() which expresses the probability for a cylinder to adopt a solid angle  on the 2D unit sphere. The shape of the ODF allows us to distinguish between isotropic order, where f = 1/4π , and nematic order where f is some peaked function. Particle interactions are incorporated on the second-virial level via a spatial integral over the Mayer function (r; 1 , 2 ) = e−βU (r;1 ,2 ) − 1,

(2)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

024901-3

H. H. Wensink and E. Trizac

J. Chem. Phys. 140, 024901 (2014)

which depends on the pair potential U between two cylinders with centre-of-mass distance r = r1 − r2 . In our model, we shall assume each particle to consist of a cylindrical hard core (HC) with diameter D and height L supplemented with a soft interaction potential Us describing (effective) long-ranged interaction with neighboring particles. These soft interactions can either be repulsive or attractive and may originate from effective interparticle forces between the colloids under the influence of depletion agents,17 polymers end-grafted onto the colloid surface,65 or electrostatics.14 The corresponding potential is unlikely to be a simple radially symmetric function but rather emerges from a particular spatial arrangement of interaction sites located on the cylinder surface. In the latter case, the soft potential is given by a summation over site-site interactions which are assumed to have a radially symmetric form u(r),66  Us (r; 1 , 2 ) = u(|r + sl (1 ) − sm (2 )|), (3)

r∈vexcl

whereas a0 represents an integration over the entire spatial volume V ,  drUs (r; 1 , 2 ) a0 = V

=

where sl denotes the distance vector between site l located on the surface of cylinder 1 and the centre-of-mass r1 . The total pair potential thus reads  ∞ if hard cores overlap (4) U (r; 1 , 2 ) = Us (r; 1 , 2 ) otherwise. For hard cylinders (Us = 0), the spatial integral over the Mayer function yields the excluded volume between two cylinders at fixed orientations. In the limit of maximal cylinder anisotropy, one obtains2  vexcl (γ ) = − drH C (r; 1 , 2 ) = v0 | sin γ |, (5) V

with v0 = 2L2 D for needles (L/D → ∞) and v0 = π D 3 /2 for disks (L/D↓0). γ (1 , 2 ) denotes the enclosed angle between the normal vectors of two cylinders. The total free energy of the fluid can be compactly written as

(6)

The spatial integral in the final term runs over the space complementary to the finite excluded volume manifold formed by the hard cores of two cylinders at fixed orientations. The last term can be interpreted as an effective excluded volume but a direct calculation of this quantity poses some serious technical difficulties.67 A more tractable expression can be obtained by adopting a mean-field form which can be obtained by taking the limit βUs 1 in the second-virial term. Equation (6) can then be recast into a form resembling a generalized van der Waals free energy ρ βF ∼ ln Vρ + ln 4πf () + vexcl (γ ) N 2 βρ (a0 − a1 (1 , 2 )) , (7) + 2 where the contributions a0 and a1 can be identified as van der Waals constants emerging from spatial averages of the soft

 l,m

dru(|r + sl (1 ) − sm (2 )|).

(9)

V

Introducing a linear coordinate transformation y → r + sl (1 ) − sm (2 ) (with Jacobian unity) yields a trivial constant  ∞  dyu(|y|) = 4π drr 2 u(r) = cst, (10) a0 = l,m

l,m

ρ βF ∼ ln Vρ + ln 4πf () + vexcl (γ ) N 2   ρ −βUs (r;1 ,2 ) dr(1 − e ) . + 2 r∈v / excl

potential. The non-trivial one, a1 , runs over the excluded volume manifold of the cylinders  drUs (r; 1 , 2 ), (8) a1 (1 , 2 ) =

V

0

independent of the mutual cylinder orientation. In arriving at Eq. (10), we have tacitly assumed that the spatial integral over the soft part of the pair potential is bounded. For this to be true, the site potential must be less singular than 1/r3 such that its 3D Fourier transform (FT) exists,  ∞ sin qr ˆ u(r). (11) u(q) = 4π drr 2 qr 0 Steep repulsive potentials such as the repulsive Coulomb (u ∼ r−1 ) or the attractive van der Waals dispersion potential68 (u ∼ −r−6 ) do not qualify and our treatment is therefore limited to cases such as the screened-Coulomb (Yukawa) potential14 or various bounded potentials such as Gaussian,69, 70 square-well,71, 72 or linear ramp potentials which routinely arise from free-volume type theories for depletion interactions17 or as effective potentials for end-grafted polymers.65 We remark that the free energy Eq. (7) represents a hybrid between the second-virial approach, which is valid at low particle densities, and the mean-field approximation, accurate at high particle density. For charged cylinders, it will be shown that the theory represents a simplified alternative to a more formal variational hard-core PB theory for anisometric colloids developed by Lue and co-workers.73, 74 We shall now proceed with analyzing the non-trivial vander-Waals contribution Eq. (8). In view of the existing FT, it is expedient to recast the spatial integral in Eq. (7) in reciprocal space. The analysis is further facilitated by using the linear transform introduced right after Eq. (9). After some rearranging the angle-dependent van der Waals term Eq. (8) can be factorized in Fourier space in the following way:  1 ˆ dqu(q)W (q; 1 )W (−q; 2 ) a1 (1 , 2 ) = (2π )3 ×vˆexcl (q; 1 , 2 ),

(12)

in terms of the FT of the excluded volume manifold of two cylinders (calculated in the Appendix),  dreiq·r vˆexcl (q; 1 , 2 ) = r∈vexcl

= v0 | sin γ |F(q; 1 , 2 ),

(13)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

024901-4

H. H. Wensink and E. Trizac

J. Chem. Phys. 140, 024901 (2014)

where the expressions for F are given explicitly in the Appendix. The contribution W pertains to a FT of the spatial resolution of the interaction sites according to  eiq·sl (α ) , α = 1, 2 (14) W (q; α ) = l

which may be interpreted as a form factor reflecting the internal structure of the interaction sites on the particle surface. The simplest case, a point segment located at the centre-ofmass thus corresponds to s1 = s2 = 0 so that W = 1. More realistic configurations shall be considered in Sec. III. Next, the equilibrium form of the ODF is obtained by a formal minimization of Eq. (7),  δ βF − λ1 = 0, (15) δf N where the Lagrange parameter λ ensures the ODF to be normalized on the unit sphere. The associated self-consistency equation for the ODF reads f (1 ) = Z −1 exp[−ρ(vexcl (1 , 2 ) − βa1 (1 , 2 ))2 ], (16) with normalization constant Z = exp[·]1 . It is easy to see that the isotropic solution f = cst, i.e., all orientations being equally probable, is a trivial solution of the stationarity condition. Beyond a critical particle density, non-trivial nematic solutions will appear which can be obtained by numerically solving Eq. (16).75 Once the equilibrium ODF is established for a given density phase equilibria between isotropic and nematic states can be investigated by equating the pressure P and chemical potential μ in both states. These are obtained by standard thermodynamic derivatives of the free energy Eq. (7), ρ2 vexcl (1 , 2 ) + βa0 − βa1 (1 , 2 ) , 2 βμ = ln ρV + ln 4πf ()

βP = ρ +

+ρ vexcl (1 , 2 ) + βa0 − βa1 (1 , 2 ) .

(17)

The thermodynamic properties of the isotropic-nematic transition can be calculated by numerically solving these coexistence equation in combination with Eq. (16), the stationarity condition for the ODF. Collective orientation order of cylinders with orientation unit vector uˆ order can be probed by introducing a common nematic director nˆ and defining nematic order parameters such as ˆ Sn = Pn (uˆ · n),

(18)

where Pn represents a nth-order Legendre polynomial (e.g., P2 (x) = (3x 2 − 1)/2). Odd contributions of Sn are strictly zero for non-polar phases and S2 is routinely used to discriminate isotropic order (S2 = 0) from uniaxial nematic order S2 = 0. III. GENERALIZED SCREENED-COULOMB POTENTIAL FOR CYLINDERS

In this section, we shall consider a simple model for charged anisotropic colloidal particles. Let us consider two disk-shaped macro-ions with total surface charge Z in an

electrolyte solution with ionic strength determined by the counter ions and additional co- and counter ions due to added salt. Formally, the electrostatic potential around the charged surface of a macro-ion in an ionic solution with a given ionic strength can be obtained from the nonlinear PoissonBoltzmann (PB) equation.14 This theory neglects any correlations between micro-ions and assumes the solvent to be treatable as a continuous medium with a given dielectric constant. In the Debye-Hückel approximation, valid if the electrostatic potential at the macro-ion surface is smaller than the thermal energy, the PB equation can be linearized and the electrostatic interaction between two point macro-ions with equal charge ±Ze in a dielectric solvent with relative permittivity εr is given by the screened-Coulomb or Yukawa form βu0 (r) = Z 2 λB

e−κr , r

(19)

with ε0 the dielectric permittivity in vacuum, r the distance between the macro-ions, λB = βe2 /4π ε0 εr the Bjerrum length (λB = 0.7 nm for water at T = 298 K), and κ −1 the Debye screening length which measures the extent of the electric double layer. In the limit of strong electrostatic screening, the screening factor is proportional to κ = (8π λB ρ 0 )1/2 with ρ 0 the concentration of added 1:1 electrolyte. In general, for highly charged colloids nonlinear effects of the PB equation can be accounted for by invoking a cell approximation76 which assumes a fully crystalline structure where each particle is compartmentalized in Wigner-Seitz cells or a so-called Jellium model77 where a tagged particle is exposed to a structureless background made up by its neighboring particles. Both methods allow for a solution of the full nonlinear PB equation for an isolated colloidal with the effect of the surrounding charged particles subsumed into a suitable boundary condition. This procedure yields so-called effective values for the charge Zeff < Z and Debye screening constant κeff which can be used to achieve accurate predictions for the thermodynamic properties (e.g., osmotic pressure) of fluids of highly charged spheres.78 We will briefly touch upon these effective parameters in Sec. III C. We reiterate that we focus here on the high-salt regime where use of the linearized form Eq. (19) combined with effective electrostatic parameters is deemed appropriate. The low-screening regime requires a lot more care due to the fact that the effective interaction becomes inherently dependent on the macroion density. As a consequence, the free energy contains non-trivial volume terms which may have important implications for the fluid phase behaviour.79–81 The FT of the Yukawa potential is given by a simple Lorentzian ˆ u(q) = Z 2 λB

4π . q2 + κ2

(20)

The spatial average over Eq. (19) yields for a0 βa0 = 4π Z 2 λB κ −2 .

(21)

We may generalize the screened-Coulomb potential for a cylindrical object by imposing the total effective electrostatic potential given by a sum over n identical Yukawa sites located on the cylinder surface. As per Eq. (3), the pair potential is

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

024901-5

H. H. Wensink and E. Trizac

J. Chem. Phys. 140, 024901 (2014)

given by

main axis of the rod pair via

βUs =

Z 2 λB  exp[−κ|r + si (1 ) − sj (2 )| . n2 i,j D. This opens up ways to designing optimized schemes that combine an effective particle shape with an appropriately rescaled aligning background potential capturing the far-field electrostatics at high particle density. These ideas have been pursued in detail in Refs. 73, 74, and 67 and shall not be further discussed here. Let us now turn to the case of charged disks. The isotropic-nematic phase diagram emerging from the Onsagervan der Waals theory for the various charge patterns depicted in Fig. 1 is shown in Fig. 4. Similar to the case of rods we observe a marked weakening of nematic order and a narrowing of the biphasic gap. The overall shape of the binodals does not depend too sensitively on the amplitude provided that Z 2 λB /D ∼ O(10) at most. Similar to what is observed for rods the isotropic-nematic ceases to exist below a critical ionic strength. This effect is most noticeable for disks with a continuous charge distribution along the face or rim. For disks with a discrete hexagonal charge patterns, the window

J. Chem. Phys. 140, 024901 (2014)

D

024901-8

30 20

10

10 − 0.15

0.6

0.65

0.7

0.75

− 0.1

0.8

− 0.05

0

S4

0.85

0.05

0.1

0.9

S2 FIG. 4. (a) Isotropic-nematic binodals for charged disks with charge distributions corresponding to the patterns indicated in Fig. 1. Plotted is the particle concentration c = ρD3 versus the ionic strength κD. The Yukawa amplitude is Z2 λB /D = 10. (b) Orientational order parameters of the nematic phase at coexistence.

of stable nematic order is somewhat larger in terms of ionic strength. The curvature of the binodals point to a reentrant phase separation phenomenon similar to the case of rods in Fig. 3. For highly charged disks, the suppression of nematic order is even more drastic and is borne out from the secondvirial free energy Eq. (35) using the orientation-dependent Yukawa potential Eq. (33). No stable isotropic-nematic was found in the experimentally relevant range of disk diameters 35λB < D < 200λB and densities. The lack of nematic stability can be inferred from Fig. 5 illustrating the effective excluded volume −β 1 of a charged disk. Although the volume depends strongly on the ionic strength its angular variation

FIG. 5. Effective excluded volume −β 1 (γ ) between highly charged disks with diameter D = 35λB interacting via the orientation-dependent Yukawa potential Eq. (33). The black solid line indicates the bare excluded volume of hard disks.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

024901-9

H. H. Wensink and E. Trizac

J. Chem. Phys. 140, 024901 (2014)

remains very weak throughout. The effective shape of a highly charged colloidal disk resembles that of a slightly deformed spherical object whose anisometry is insufficient to generate a thermodynamically stable orientational disorder-order transition. Like for the rod case, improved schemes could be envisaged for the strong-coupling regime by combining an effective disk shape with a suitably chosen amplitude for the mean-field aligning potential a1 (γ ). V. STABILITY OF LIQUID CRYSTAL PHASES WITH POSITIONAL ORDER

Possible phase transitions to spatially inhomogeneous states with smectic or columnar order can be investigated by recasting the mean-field Onsager into a functional form depending on the one-body density field ρ(r, ).92 Within the framework of classical DFT, the free energy functional needs to be minimized with respect to ρ to yield the unique equilibrium density profile for a given chemical potential, temperature, and external potential.93 In this work, we shall adopt a simple stability analysis94, 95 by assuming a weak periodic density modulation with wave-vector k and amplitude ε, ρ(r, ) = ρ0 f0 () + εf ∗ () cos(k · r),

(37)

superimposed onto the one-body density ρ 0 (r, ) = ρ 0 f0 () of the spatially homogeneous phase. Beyond a particular value of the bulk density such a periodic density perturbation will lead to a reduction of the free energy and the homogeneous bulk phase will become marginally unstable. The so-called bifurcation point can be found by inserting Eq. (37) into the density functional and Taylor-expanding up to second order in ε. The resulting bifurcation condition is represented by a linear eigenvalue equation,96  ∗ ˆ f (1 ) = ρ0 f0 (1 ) d2 f ∗ (2 )(k; 1 , 2 ), (38) in terms of the cosine transformed Mayer function  ˆ (k; 1 , 2 ) = dr(r; 1 , 2 ) cos(k · r).

(39)

The eigenvector f*() probes the angular distribution in the new phase and reflects the intrinsic coupling between positional and orientational order. A bifurcation to the positionally modulated state occurs at the wave vector k that generates the smallest eigenvalue ρ 0 > 0 of Eq. (38). If we neglect the translation-rotation coupling, then f*() = f0 () and the bifurcation condition takes the form of a divergence of the static structure factor S(k), ˆ 1 , 2 )) = 0. S(k)−1 = (1 − ρ0 (k;

(40)

By applying the van der Waals approximation outlined in ˆ can be expressed as a sum of the hard-core conSec. II,  tribution and a part that encodes the effect of the soft potential. Eliminating the angular dependency of the excluded volume vˆexcl and form factor W for notational brevity, one arrives at the following expression for the Mayer kernel in Fourier

space: ˆ ˆ (k) = −vˆexcl (k) − u(k)  1 ˆ + dqu(q)W (q)W (−q)vˆexcl (k − q). (2π )3 (41) The Fourier integral presents a non-trivial mode-coupling term that convolutes the imposed density wave with the modes describing the distance-dependence of the soft interactions. The solution of Eq. (40) (or Eq. (38)) for particles with full orientational degrees of freedom poses a substantial technical task and we shall simplify matters by considering the more tractable case of parallel cylinders. Let us equate the particle frame to the lab frame {ˆx, yˆ , zˆ } with the cylinder normals pointing along zˆ . The excluded volume of two parallel cylinders is again a cylinder with volume 2π LD2 . In Fourier space, the excluded volume takes the following form:  J1 ( (Dq · xˆ )2 + (Dq · yˆ )2 ) . vˆexcl (q) = 2π LD 2 j0 (Lq · zˆ ) 1  2 + (Dq · y 2 ˆ ˆ ) (Dq · x ) 2 (42) Due to the parallel orientation it is no longer possible to take the limit of infinite particle anisometry since the excluded volume vanishes in both limits (similar to setting γ = 0 in Eq. (5)). Therefore, we shall consider the case D/L x (rods) and L/D x (disks) with x a small but finite number and use the volume fraction φ = (π /4)LD2 ρ 0 as a convenient measure for the particle concentration. We may probe instabilities pertaining to smectic order by identifying k = kS {0, 0, 1}, a one-dimensional periodic modulation along the nematic director. Hexagonal columnar order can be parametrized by a linear superposition of three modu√ 3 1 lations with wavevectors k1 = kC {0, 1, 0}, k2 = kC { 2 , 2 , 0}, √

and k3 = kC {− 23 , 12 , 0} describing a two-dimensional triangular lattice perpendicular to the director. The results in Fig. 6 reveal a marked stabilization of columnar at the expense of smectic order for rodlike cylinders at low ionic strength. This outcome is in accordance with previous numerical results for Yukawa rods in a strong external aligning field.54 Needless to say that the transition values are merely qualitative and that the volume fractions can be brought down to more realistic values, for instance, by using an effective second-virial theory based on a resummation of higher virial coefficient (e.g., using Parsons’ theory58 ). For hard parallel cylinders, the nematic-smectic always pre-empts the nematic-columnar one irrespective of the aspect ratio x. This implies that the parallel approximation fares rather badly for hard discotic systems which are known to form columnar phases only.9, 97 Nevertheless, some general trends can be gleaned from Fig. 6(b) such as an apparently stabilization of smectic order for uniformly charged disks in the low screening regime. The prevalence of smectic order has been recently reported in weakly screened discotic systems.98 As for the other charge patterns, the observation from Fig. 6 that both smectic and columnar-type order are destabilized upon reducing the screening could hint at more complicated instability mechanisms prevailing in the low screening region, such as those pertaining to crystalline order where both longitudinal

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

024901-10

H. H. Wensink and E. Trizac

J. Chem. Phys. 140, 024901 (2014)

where α is a parameter which describes how localized the particles are around each lattice site. If we assume the proportion of lattice defects to be negligible, each lattice site should contain only one particle as reflected in the normalization of Eq. (44). The excess free energy of the system can be expressed in terms of the following Fourier integral:  kB T  1 ˆ 2 (k; ˆ ,  ), dkeik·Rij G(k) Fex = − 2 i =j (2π )3 (45) in terms of the FT of the orientation-dependent Mayer kernel ˆ Eq. (41), and the Gaussian weight G(k) = exp(−k 2 /4α). In general, the radially symmetric form Eq. (44) is justified only if particles are strongly localized around their lattice points (α  1) so that the density peaks are not affected by the symmetry of the underlying lattice. The total free energy is obtained by combining the excess free energy with the ideal free energy associated with the Gaussian parameterization  2   V α 5 3 ln − + ln f () . (46) Fid = N kB T 2 π 2

FIG. 6. (a) Variation of the nematic-smectic (NS) and nematic-columnar (NC) bifurcation density with ionic strength κD for parallel charged rods with aspect ratio L/D = 50 and (b) parallel disks with surface charge patterns indicated in Fig. 1 (D/L = 10, Z2 λB /D = 10). Solid curves indicate NS bifurcations, dotted curves NC instabilities.

and transverse density modulations compete with spatial inhomogeneities in the director field.47 We wish to underline that the approach outlined above is generic in that it provides a simple route to gauge the effect of soft interactions on the stability of positionally ordered liquid crystals. It can be applied to a vast range of model systems with various segment potentials (provided integrable) and form factors. Instabilities from nematic to other liquid crystals symmetries or three-dimensional crystals (e.g., fcc or bcc) can easily be included by adapting the k vectors to the desired Bravais lattice. In order to describe fully crystalline states, we may exploit the fact that particles are strongly localized around their lattice site to construct an appropriate density functional representation for the excess Helmholtz free energy. In the following, we shall briefly sketch the approach outlined in Refs. 99 and 100. The central assumption is that the density profile of the solid consists of Gaussian peaks centred on a predefined lattice vector {Ri } factorized with the orientational probability (ODF) f(). If we assume a spatially homogeneous director field, the one-body density can be written as ρ(r, ) = f ()

N 

G(r − Ri ),

(43)

exp[−α(r − Ri )2 ],

(44)

i=1

with G(r − Ri ) =

 α 3/2 π

Next, the free energy must be minimized with respect to the localization parameter α, the set of relevant lattice constants corresponding to the imposed lattice symmetry101 and f(). This simple variational scheme allows one to compare the stability of various crystal symmetries as a function of density and interaction range and strength. In addition, due to the translation-orientation coupling via f both aligned and rotationally disordered, plastic crystal states can be included. Phase transitions between fluid and crystal phases can be probed by equating the pressure and chemical potential emerging from the Gaussian free energy with those of the fluid phases Eq. (17). VI. CONCLUDING REMARKS

We have proposed a generalized Onsager theory for strongly non-spherical colloidal particles with an intrinsic patchiness in the interaction potential. The theory supplements the second-virial reference free energy for the hardcore interaction with a first-order perturbative (van der Waals) term which captures the directional soft interactions between the rods or the disks. As such, the theory interpolates between the low density regime, where the second-virial approximation holds, and the high density regime where the mean-field approach is accurate. We have aimed at formulating a generic framework that should be applicable to a wide range of particle shapes, ranging from elongated rods to flat, sheet-like disks with an arbitrary spatial organization of interactions sites distributed along the colloid surface. By recasting the mean-field contribution in terms of a Fourier series, the excess free energy naturally factorizes into three main contributions; the site-site interaction potential, the shape of the colloidal hard-core, and a form factor associated with the spatial arrangement of the interaction site residing on each particle. As a test case, we have applied our theory to investigate orientation disorder-order transitions in fluids of charged rods and disks with a uniform, localized, or discretized surface charge pattern. The results for the isotropic-nematic phase

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 70.119.188.32 On: Tue, 25 Mar 2014 17:04:55

024901-11

H. H. Wensink and E. Trizac

diagram and the instability analysis of transverse and longitudinal freezing of a nematic fluid in the high-density regime reveal a picture that is consistent with results from experiment and particle simulation. Evidence for the trends predicted by our theory can be drawn from various experimental observations that we will summarily discuss next. A reduction of the biphasic gap as well as indications of a non-monotonic dependence of the isotropic- (cholesteric) nematic phase boundaries upon decreasing ionic strength have been reported in systems of stiff fd rods.102 A similar reduction of the phase gap was reported for rigid tobacco mosaic virus (TMV) rods103 and high aspectratio zirconium phosphate platelets.104 The marked destabilization and indeed complete absence of nematic order for charged platelets at low ionic strength (Fig. 4) is in line with results from recent simulation of particles with a hard-core Yukawa potential.47 In addition, charged gibbsite platelets cease to form stable nematic phases at low ionic strength in favor of columnar order.105 As for the high-density regime discussed in Sec. V, a crossover from smectic order to a more intricate ordered state upon decreasing ionic strength has been observed in concentrated systems of TMV rods.106 However, it is not fully clear whether these structures are really columnar or represent three-dimensional crystalline order. A much more convincing account of the preference of hexagonal columnar over smectic order of elongated charged rods (see Fig. 6(a)) has been reported for semiflexible fd virus rods.107 Conversely, the tendency of charged platelets to self-assemble into layered structures, hinted at by the results in Fig. 6(b), was highlighted in recent studies of charged niobate nanosheets108 and gibbsite platelets at low-screening solvent conditions.98 Recent simulation work on discotic systems with explicit point charges demonstrates that similar layered, smectic phases may be formed by oppositely charged oblate mesogens.42 These observations lend credence to our theory as a practical tool to assess the influence of soft patchy interactions on the liquid crystal phase diagram of non-isometric colloids. Although the focus of this study is on the liquid crystal fluid phases that emerge at relatively low particle density, the stability of spatially ordered liquid crystals at higher particle concentration can also be scrutinized using a simple bifurcation analysis while fully crystalline phases can be expediently accounted for using a Gaussian parameterization for the one-body density often used in density functional theories of freezing. We remark that the present theory is amenable to various extensions towards more complicated systems. Colloidal dispersions composed of non-spherical particles are rarely monodisperse but are often characterized by a continuous spread in particle sizes. The polydisperse nature of the colloid shape and/or the amplitude of the soft interactions can be incorporated in a straightforward manner.109, 110 Bio-colloids such as stiff viral rods111 and DNA are commonly characterized by an intrinsic helical patchiness which has profound implications on the mesostructure in bulk and confinement.112 The present theory could be extended to relate the mesoscopic chirality of twisted nematics to the intrinsic helical form factor of the colloid.113

J. Chem. Phys. 140, 024901 (2014)

Last but not least, similar to systems of spherical subunits,114, 115 more accurate reference free energies could be employed which should give a more reliable account of correlations in systems of less anisometric colloids (dumbbells, thick platelets, polyhedra) which routinely form highly ordered (liquid) crystals at high particle volume fraction.116, 117 APPENDIX: EXCLUDED VOLUME OF STRONGLY ANISOMETRIC CYLINDERS

In this appendix, we derive expressions for the FT of the excluded volume manifold of two infinitely slender rods and disks, featured in Eq. (12) of the main text. The excluded volume of two hard cylinders at fixed angle γ is a parallelepiped which can be parameterized by switching from the laboratory frame to a particle frame spanned by the normal orientational unit vectors uˆ α of the cylinder pair. Let us define the additional unit vectors vˆ | sin γ | = uˆ 1 × uˆ 2 , ˆ α = uˆ α × vˆ w

(α = 1, 2),

(A1)

ˆ α } are two orthonormal basis sets in 3D. The so that {uˆ α , vˆ , w centre-of-mass distance vector can be uniquely decomposed in terms of these basis vectors ˆ α )w ˆ α. r = (r · uˆ α )uˆ α + (r · vˆ )ˆv + (r · w

(A2)

The leading order contribution to the excluded-volume body is of O(L2 D) and stems from the overlap of the cylindrical parts of the spherocylinders. This resulting parallelepiped can be parameterized as follows: rcc =

L L t1 uˆ 1 + t2 uˆ 2 + Dt3 vˆ , 2 2

(A3)

with −1 ≤ ti ≤ 1 for i = 1, 2, 3. The Jacobian associated with the coordinate transformation is Jcc = 14 L2 D| sin γ |. The FT of the parallelepiped is thus given by  vˆexcl (1 , 2 ) = drcc eiq·rcc = Jcc

 i