Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Coarse-grained treatment of the self-assembly of colloids suspended in a nematic host phase

Sergej Püschel-Schlotthauer a, Tillmann Stieger a, Michael Melle a, Marco G. Mazza b and Martin Schoen ac
aStranski-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
bMax-Planck-Institut für Dynamik und Selbstorganisation, Am Faßberg 17, 37077 Göttingen, Germany
cDepartment of Chemical and Biomolecular Engineering, North Carolina State University, Engineering Building I, Box 7905, 911 Partners Way, Raleigh, NC 27695, USA

Received 27th July 2015 , Accepted 9th October 2015

First published on 13th October 2015

The complex interplay of molecular scale effects, nonlinearities in the orientational field and long-range elastic forces makes liquid-crystal physics very challenging. A consistent way to extract information from the microscopic, molecular scale up to the meso- and macroscopic scale is still missing. Here, we develop a hybrid procedure that bridges this gap by combining extensive Monte Carlo (MC) simulations, a local Landau–de Gennes theory, classical density functional theory, and finite-size scaling theory. As a test case to demonstrate the power and validity of our novel approach we study the effective interaction among colloids with Boojum defect topology immersed in a nematic liquid crystal. In particular, at sufficiently small separations colloids attract each other if the angle between their center-of-mass distance vector and the far-field nematic director is about 30°. Using the effective potential in coarse-grained two-dimensional MC simulations we show that self-assembled structures formed by the colloids are in excellent agreement with experimental data.

1 Introduction

Liquid crystals are fluids made of molecules that lack spherical symmetry. Instead, their molecules contain elongated, rigid cores that form nematic liquid crystals, or disk-like cores that produce discotic liquid crystals or even more complex shapes. This simple property of anisotropy produces a myriad consequences for the optical, mechanical and thermodynamic properties of liquid crystals. For example, as the temperature is decreased they undergo a series of phase transitions where the symmetry of their state is spontaneously broken. Starting from a high temperature isotropic fluid, where all orientations are equivalent, to a nematic state where orientational order is broken, to a smectic phase, where positional order is broken in one dimension.

The molecular anisotropy produces effects on a macroscopic scale. In the nematic state, all molecules tend to align in the same direction, called the fluid's director.1 But any localized deviation of molecular orientations from the director will cause a restoring, elastic force. However, a single global orientation cannot be satisfied for all boundary conditions. A colloid immersed in the nematic fluid causes a specific preferential alignment of the liquid crystal molecules on its surface, which is termed anchoring. A spherical colloid has a symmetry incompatible with a global nematic director. Thus, defects in the orientational field will arise. These topological defects are points or lines that minimize the free energy of the liquid crystal subject to complex boundary conditions. Depending on details of the fluid and the anchoring of its molecules at the colloid, the orientational field can be of such dazzling complexity that experts are just beginning to unravel its structural details.2

If multiple colloids are immersed in a nematic liquid crystal, the distortions of the local director field also give rise to effective interactions between the colloids mediated by the nematic host fluid.3 These interactions may be used to self-assemble the colloids into supramolecular entities in a controlled (i.e., directed) manner. In this way ordered assemblies of colloids of an enormously rich structure with a great variety of symmetries may be built that would not otherwise exist without the ordered nature of the host phase.4,5 These self-assembled colloidal structures are also of practical importance, as they can be used to produce photonic crystals.6,7

The forces guiding colloidal self-assembly result from a complex interplay of molecular scale ordering, mesoscopic elastic interactions, and large scale topological arguments. While the framework of the Landau–de Gennes theory is appropriate for mesoscale effects, the elastic Frank–Oseen free energy is appropriate for long-range interactions.8 However, there is a gap between the microscopic, molecular information and the coarse-grained approaches used at the meso and macroscopic scale. No consistent method exists to transfer physical information from the molecular scale up to the scale where elastic forces and nonlinearities in the nematic interactions occur. Here, we pursue the goal of bridging this gap. The task is not idle, as there are physical situations that have so far eluded a precise explanation at the molecular level. As a test case to demonstrate power and validity of a novel hybrid approach developed in this work we take on the problem of colloidal self-assembly. Poulin and Weitz9 found experimentally that in a nematic phase the colloidal center-to-center distance vector r12 forms a “magic” angle of θ ≈ 30° with the global director [n with combining circumflex]0 if the liquid crystal molecules (mesogens) favor a local planar anchoring at the colloid's surface.10

Near an isolated colloid with planar anchoring the mesogens will produce a Boojum defect. While for such a single topological defect long-range approximation analogous to electrostatics can be sufficient, the complex interaction of multiple Boojum defects requires a more careful analysis. In previous theoretical attempts a much larger angle of about 50° is usually found.9,11 This number is based upon calculations where one employs the electrostatic analog of the Boojum defect topology.9 In fact, as stated explicitly by Poulin and Weitz “This theoretical value is different from the experimentally observed value for θsince the theory is a long-range description that does not account for short-range effects”.9

Theoretical results similar to the experimental ones by Poulin and Weitz9 were recently found by Tasinkevych et al. through a direct minimization of a Landau–de Gennes free-energy functional.8 These authors demonstrate that the radial component of the elastic force has an attractive minimum around θ ≈ 30° for certain r12 = |r12|; this minimum shifts to larger θ as r12 increases. However, to date and to the best of our knowledge no molecular-scale work exists supporting the result of Tasinkevych et al.8 or the experimental findings by Poulin and Weitz.9

Another motivation for our work is the more recent experimental observation that between a pair of colloids in a nematic host repulsive and attractive forces act depending on θ.11 For example, at θ ≈ 30° the colloids attract each other whereas at θ = 0° and 90° repulsion between the colloids is observed.

To study these effects starting from a molecular-level based description we employ a combination of Monte Carlo (MC) simulations in the isothermal–isobaric and canonical ensembles, two-dimensional (2D), coarse-grained MC simulations in the canonical ensemble, classical density functional theory (DFT), concepts of finite-size scaling (FSS), and Landau–de Gennes (LdG) theory to investigate the effective interaction between a pair of spherical, chemically homogeneous colloids mediated by a nematic host phase.

The remainder of this manuscript is organized as follows. We begin in Section 2 with an introduction of our model system and its various constituents. Relevant theoretical concepts are introduced in Section 3. In Section 4 we present results of this study which we summarize and discuss in the closing Section 5 of this manuscript.

2 Model

In this work we consider the self-assembly of colloidal particles in two and three dimensions. The colloids are immersed in a nematic liquid-crystalline host fluid where an external field is invoked to control the global director [n with combining circumflex]0. The following three sections are devoted to introducing the various constituents of our model and to specifying their interactions with one another.

2.1 Host phase

The liquid crystal host phase is composed of N mesogens interacting with each other in a pairwise additive fashion. The interaction potential can be cast as
φmm(rij, ωi, ωj) = φiso(rij) + φanis(rij, ωi, ωj)(1)
where rij = rirj is the distance vector connecting the centers of mass of a mesogenic pair located at points ri and rj, respectively, and rij = |rij|. According to eqn (1) the full interaction potential is split into an isotropic (φiso) and an anisotropic part (φanis) where the latter depends on the orientations ωi and ωj of the mesogenic pair. In fact, ωi,j = (θi,j, ϕi,j) where θi,j and ϕi,j are Euler angles implying that the mesogens have uniaxial symmetry.

For the isotropic contribution we adopt the well-known Lennard-Jones potential

image file: c5sm01860a-t1.tif(2)
where εmm is the depth of the attractive well, σ is the van der Waals diameter of the isotropic core, and φrep and φatt represent repulsive and attractive contributions to φiso, respectively.

To derive an expression for the anisotropic contribution in eqn (1) we follow Giura and Schoen.12 From a lengthy but relatively straightforward derivation these authors show that

φanis(rij, ωi, ωj) = φatt(rij)Ψ([r with combining circumflex]ij, ωi, ωj)(3)
In eqn (3) the anisotropy function is given by
Ψ([r with combining circumflex]ij, ωi, ωj) = 5ε1P2[û(ωiû(ωj)] + 5ε2{P2[û(ωi[r with combining circumflex]ij] + P2[û(ωj[r with combining circumflex]ij]}(4)
where û(ωi,j) and [r with combining circumflex]ij = rij/rij are unit vectors specifying the orientation of mesogens i and j and the orientation of the center-of-mass distance vector, respectively, both in a space-fixed frame of reference. In addition, P2(x) = ½(3x2 − 1) is the second Legendre polynomial and the (dimensionless) anisotropy parameters 2ε1 = −ε2 = −0.08 are fixed throughout this work. Notice also that the first term on the right-hand side of eqn (4) matches the orientation dependence of interactions in the Maier–Saupe model13,14 whereas the next two terms account for the anisotropy of mesogen–mesogen attractions with enhanced sophistication.

2.2 External field

For a suitable choice of thermodynamic state parameters the host phase introduced in Section 2.1 exhibits isotropic and nematic phases.12 Unfortunately, when a nematic phase is forming in the bulk it is impossible to determine beforehand the direction of [n with combining circumflex]0. In fact, there is an infinite number of easy axes15 with which [n with combining circumflex]0 may align in the bulk nematic phase.

To predict and control [n with combining circumflex]0 it has therefore become customary to place the liquid crystal between solid substrates with specially prepared surfaces that tend to align mesogens in their vicinity in a desired way.16 Because orientational order in a nematic liquid crystal is a long-range phenomenon, substrate-induced alignment of mesogens allows one to control [n with combining circumflex]0 on a length scale exceeding by far a molecular one. A host of different techniques to achieve a particular substrate anchoring of mesogens experimentally has been known for quite some time.17

In this work we take the substrates to be planar and structureless such that their surfaces are separated by a distance sz along the z-axis. Thus, the interaction between a mesogen and the substrates can be cast as

image file: c5sm01860a-t2.tif(5)
where Δzi = zi ± sz/2 is the distance of mesogen i from the lower (+) and upper (−) substrate plane, respectively. Hence, φms may be viewed as a local, orientation dependent external field, the strength of which is controlled by εms. Throughout this work we take εms/εmm = 5.00.

The value of εms is chosen for two reasons. First, it guarantees a sufficiently strong alignment of mesogens with the surface so that fluctuation of [n with combining circumflex]0 over the course of the simulations are negligible. Second, εms is still small enough to prevent from freezing those portions of the host phase that are located in the vicinity of the substrates.

The orientation dependence of the external field in our model enters through the anchoring function

g(ωi) = [û(ωiêx]2(6)
where êx is a unit vector parallel to the x-axis. Hence, the anchoring function discriminates energetically between a desired orientation of mesogens parallel with this axis and less desired ones for which |û(ωiêx| < 1. In other words, g(ωi) may be viewed as a mathematical “device” mimicking aligning substrates in experimental setups. On account of its definition in eqn (6), g(ωi) allows one to more or less pin [n with combining circumflex]0 to the x-axis on average where |[n with combining circumflex]0·êx| ≃ 1 due to thermal fluctuations.

2.3 Colloidal particles

In addition, colloidal particles are immersed into the nematic liquid crystal. These colloids are spherical in shape, where r0 = 3.00σ is their hard-core radius, and have a chemically homogeneous surface. Similar to the planar substrates we envision the surfaces of the colloids to have been treated such that mesogens have a specific orientation with respect to the local surface normal of a colloid. Following earlier work by some of us18 we take the mesogen–colloid interaction potential to be given by
image file: c5sm01860a-t3.tif(7)
where rcij = |rircj| is the distance between the center of mass of mesogen i at ri and that of colloid j at rcj. Hence, rcijr0 is the distance of the center of mass of mesogen i from the surface of colloid j. The strength of the mesogen–colloid interaction is controlled by εmc which we maintain constant so that εmc/εmm = 3.50. Again, this value has been selected on the basis of the same rationale given in Section 2.2.

In eqn (7), η is the inverse Debye screening length and parameters a1 and a2 have been introduced to guarantee that the minimum of φmc remains at a distance r0 + σ from the colloid's center of mass and to maintain the depth of the attractive well at εmc irrespective of η.18 Throughout this work we fix ησ = 0.50.

Last but not least, we introduce another anchoring function in eqn (7) which we take to be given by

gc([r with combining circumflex]cij, ωi) = [1 − |[r with combining circumflex]cij·û(ωi)|]2(8)
where [r with combining circumflex]cij = rcij/rcij is the local surface normal of the colloid. Hence, the anchoring function introduced in eqn (8) serves to align mesogens in a degenerate,15 locally planar fashion with respect to the colloid's surface.

3 Theory

To eventually simulate the self-assembly of several colloidal particles in a coarse-grained fashion we seek to represent the nematic host phase through an effective interaction potential. Key to this approach (as in all coarse-grained treatments) is a sensible protocol to integrate out irrelevant degrees of freedom while preserving the correct physics. Here, we seek to link a molecular-level description of the nematic host to a continuum treatment.

3.1 Continuum approach

As already pointed out in Section 1, the presence of a colloid in a nematic liquid crystal causes [n with combining circumflex](r) to deviate from [n with combining circumflex]0 in certain regions near the colloid's surface. Associated with this distortion of [n with combining circumflex](r) is a local deviation between the nematic order parameter and its bulk value in the absence of any colloid. Both, the distortion of [n with combining circumflex](r) and the associated decline of S(r) cause changes in the free energy of the composite system (i.e., colloid plus nematic host). Adopting a continuum perspective a key quantity is the local alignment tensor Q(r) whose components can be cast as
image file: c5sm01860a-t4.tif(9)
where S(r) is the nematic order parameter, nα(r) is the α-component of [n with combining circumflex](r), and δαβ is the Kronecker symbol. The assumption underlying eqn (9) is that a spatial variation of the degree of nematic order of uniaxial symmetry and a deformation of the nematic director field are coupled (see also Section 4.4).

Therefore, the total change in free-energy density may readily be expressed as

Δf(r) = ΔfLdG(r) + fel(r) + fcore(10)
where the Landau–de Gennes contribution is given by
image file: c5sm01860a-t5.tif(11)
using Einstein's summation convention. Eqn (11) is essentially a Taylor expansion of the free-energy density in terms of the order parameter tensor at the isotropic-nematic (IN) phase transition. In eqn (11), A, B, and C are unknown expansion coefficients depending only on density ρ and temperature T. Assuming uniaxiality, we employ the identities QαβQβα = [/]S and QαβQβγQγα = ¾S3 This allows us to rewrite eqn (11) as
image file: c5sm01860a-t6.tif(12)
To simplify the notation we also temporarily dropped the argument r of the components of Q(r) on the righthand side of eqn (11).

In eqn (11) as well as in eqn (12)

Δf0 = ASb2 + BSb3 + CSb4(13)
is a similar LdG free-energy density of the nematic host phase obtained for the same T and ρ but in the absence of colloids and relative to the free-energy density of a corresponding isotropic fluid. In eqn (13), Sb is the (global) bulk nematic order parameter. To arrive at the expression in eqn (13) the same identities linking products of the alignment tensor to the nematic order parameter have been used as before.

Next, the elastic contribution to Δf in the one-constant approximation (see below) may be cast as

image file: c5sm01860a-t7.tif(14)
where fFO(r) is the (local) Frank–Oseen free-energy density and eqn (9) has also been used. In eqn (14), L is an elastic and K is the Frank constant. The two are related via
image file: c5sm01860a-t8.tif(15)
To apply eqn (14) some caution is advisable. This is because the expression for fFO(r) in eqn (14) is derived under the assumption that spatial variations of [n with combining circumflex](r) occur on a length scale that is large compared to a molecular one. As we shall demonstrate below this is true in our system almost everywhere except inside the core of defects. To avoid an improper calculation of fFO(r) we restrict the evaluation of eqn (14) to regions outside the defect core.19–21 The latter is defined through the inequality S(r) ≤ SIN where SIN is the nematic order parameter at the IN phase transition in the host phase (and in the absence of any colloidal particle; see below).

Because we are restricting the use of eqn (14) by “cutting out” the defect core some correction to Δf due to the core region is required. As pointed out by de las Heras et al. it is necessary to include such a correction to the change in free-energy density because of the nanoscopic size of our colloidal particle.22 This correction, represented by fcore in eqn (10), will be discussed in some detail in Section 4.3.

Last but not least, we emphasize that within the framework of the present continuum theory the standard approach is to minimize the functional23

image file: c5sm01860a-t9.tif(16)
However, this procedure has a twofold drawback. First, the minimization bears the danger that its solutions S(r) and [n with combining circumflex](r) do not necessarily correspond to the global minimum of Δ[scr F, script letter F]. This is in particular so if the structure of S(r) and [n with combining circumflex](r) in thermodynamic equilibrium is rather complex.24 Second, there is no way within the framework of the continuum approach to determine the material constants A, B, C, and K (or L). Here one usually resorts to experimental information which is available for a few liquid crystals.25

In closing, we emphasize that the expressions given in eqn (12) and (14) correspond to the same ground state. For example, in the absence of any perturbation, that is if S(r) = Sb and [n with combining circumflex](r) = [n with combining circumflex]0, ΔfLdG = fel = 0.

3.2 Molecular approach

3.2.1 Basic properties. Here we pursue an alternative approach based upon a molecular picture of the host phase and suitable for implementation in MC simulations. If carried out according to the acknowledged rules of the art, MC gives us immediate access to S(r) and [n with combining circumflex](r) for a thermodynamic equilibrium situation via suitably defined ensemble averages. One may then feed S(r) and [n with combining circumflex](r) into expressions such as eqn (12) and (14) to eventually obtain the absolute minimum of Δ[scr F, script letter F] [after an integration of Δf(r) over volume, see also Section 4.3].

Thus, to eventually compute Δ[scr F, script letter F] we begin by introducing the local alignment tensor

image file: c5sm01860a-t10.tif(17)
at the molecular level where 1 is the unit tensor and 〈…〉 denotes an ensemble average.18,26 In eqn (17), ρ(r) is the local number density. Because Q(r) can be represented by a 3 × 3 matrix it has three eigenvalues which we obtain numerically using Jacobi's method.27 We take the largest eigenvalue of Q(r) as the local nematic order parameter S(r) and the associated eigenvector as the nematic director field [n with combining circumflex](r).

However, to compute ΔfLdG(r) and fel(r) the material constants A, B, C, and K are required. Whereas this is relatively straightforward within the framework of MC simulations as far as K is concerned, one encounters serious difficulties to compute A, B, and C reliably. We address these difficulties below.

For the calculation of K we employ a method suggested by Allen and Frenkel.28,29 In their approach one considers fluctuations of Fourier components [Q with combining tilde](k) of Q(r). In the limit of |k| → 0 simple linear relationships exist from which the Frank constants K1, K2, and K3 associated with bend, twist, and splay deformations of [n with combining circumflex](r) can be estimated reliably. Stieger et al.30 have recently applied the method of Allen and Frenkel28,29 to show that for the present model of the host phase K1K2K3 = K so that the one-constant form1 of fel in eqn (14) is an excellent approximation. Under the thermodynamic conditions used here (see Section 4.1), K ≃ 1.66εmmσ−1 is obtained.

3.2.2 Classical density functional theory. To compute A, B, and C the situation is more difficult. In principle, one could obtain these constants from the order-parameter distribution P(Sb) obtained for a bulk nematic phase without colloidal particles. However, as discussed in detail by Eppenga and Frenkel26 and later by Greschek and Schoen31 it is next to impossible to determine B and C with sufficient statistical accuracy from P(Sb).

We therefore resort to a different approach based upon classical mean-field DFT. As demonstrated elsewhere32 the change in free energy–density of the bulk nematic relative to the isotropic phase can be cast as

image file: c5sm01860a-t11.tif(18)
where β = 1/kBT (kB is Boltzmann's constant), x = cos[thin space (1/6-em)]θ, and θ is the azimuthal angle if we take the z-axis of our coordinate system to coincide with [n with combining circumflex]0. Members of the set {Sl} defined on the interval [0, 1] are order parameters and {ul} are parameters that account for the contribution of anisotropic mesogen–mesogen interactions to the free energy, respectively.32

The integrand on the right-hand side of eqn (18) contains the orientation distribution function [small alpha, Greek, macron](x). It depends only on θ due to the uniaxial symmetry of the nematic phase and is normalized according to

image file: c5sm01860a-t12.tif(19)
which implies that in the isotropic phase [small alpha, Greek, macron](x) = ½. Again, because of the uniaxiality of the nematic phase we expand [small alpha, Greek, macron](x) in terms of Legendre polynomials {Pl} via
image file: c5sm01860a-t13.tif(20)
Inserting this expression into eqn (18), expanding the integrand in a Taylor series around ξ = 0 (i.e., at the IN phase transition), and retaining in this expansion terms up to [scr O, script letter O] (ξ4) allows us to recast eqn (18) as
image file: c5sm01860a-t14.tif(21)
if we limit ourselves to the leading term l = 2 in the expression for ξ and use S2 = Sb. In eqn (21), a = 2ρkB/5 and T* = −5ρu2/2kB where the latter is the temperature at which the isotropic phase becomes thermodynamically unstable. Assuming Δfn = Δf0 one easily obtains
image file: c5sm01860a-t15.tif(22a)
image file: c5sm01860a-t16.tif(22b)
by comparing terms of equal power in Sb in eqn (21) with corresponding ones in eqn (13). Hence, B < 0, C > 0, and A changes sign at T = T* as it is to be expected at the IN phase transition.1

3.2.3 Elements of finite-size scaling theory. Unfortunately, the expression for T* given in the preceding section depends on u2 which in itself depends on the level of sophistication at which pair correlations are treated within mean-field DFT. For example, at simple mean-field (SMF) level, where one completely ignores pair correlations, u2 = − 32πε1εmm/15. A more elaborate, temperature dependent form for u2 obtains at so-called modified mean-field (MMF) level [see eqn (3.7) and (3.8) of ref. 12] where one takes into account pair correlations to some extent via an orientation-dependent Mayer f-function.

Unfortunately, at SMF level T* turns out to be grossly underestimated whereas at MMF level its value is equally grossly overestimated as a previous FSS study of the IN phase transition suggests.31 This failure can be linked to the complete neglect of pair correlations at SMF level and their overestimation in liquidlike phases at MMF level.12

To improve this situation we pursue a different approach invoking concepts of FSS theory. First, within LdG theory it is a relatively simple matter to show that33

image file: c5sm01860a-t17.tif(23)
Second, with an improved estimate for the temperature TIN at the IN phase transition and coefficients for a, B, and C from DFT one can hope to obtain an improved estimate for T* from the above expression.

In FSS theory one makes explicit use of the fact that in any molecular simulation one is always confronted with systems of finite extent. Considering moments

image file: c5sm01860a-t18.tif(24)
of the order-parameter distribution the key quantity in FSS are cumulants of its various moments. Here, the second-order cumulant
image file: c5sm01860a-t19.tif(25)
is particularly useful.31 If a phase transition is discontinuous in principle (as the IN phase transition) but rounded on account of the finiteness of the system under study one anticipates pairs of cumulants for different system sizes to have a common intersection31 which scales as [script L]d where [script L] is the linear extent of a system of dimension d.34

Moreover, it has been demonstrated by Vollmayer et al.35 that the “distance” of a unique cumulant intersection from the point at which the phase transition would occur in the thermodynamic limit scales as [script L]−2d. Thus, one can expect that for sufficiently large systems it may seem that even at a discontinuous phase transition all cumulants intersect in a unique point which then for all practical purposes may be taken as the state point at which the phase transition would occur in the thermodynamic limit.

4 Results

4.1 Numerical details

Our results are based upon MC simulations in a specialized isothermal–isobaric ensemble discussed in detail elsewhere.18 In this ensemble a thermodynamic state is specified by N, T, the ratio of side lengths of the simulation cell in the x- and y-directions sx/sy, the distance sz between the solid substrates, and the transverse component P = ½(PxxPyy) of the pressure tensor P. The specialized isothermal–isobaric ensemble is employed to equilibrate the system. Production runs were carried out in the canonical ensemble using the average side lengths from the isothermal–isobaric equilibration run as fixed input values.36

We generate a Markov chain by a conventional Metropolis protocol where it is decided with equal probability whether to displace a mesogen's center of mass by a small amount or to rotate the mesogen around a randomly chosen axis. All mesogens are considered sequentially such that a MC cycle commences once each of the N mesogens has been subjected to either a displacement or rotation attempt.

Our results are typically based upon 1.5 × 105 MC cycles for equilibration followed by another 1.0 × 106 cycles during which relevant ensemble averages are taken. To save computer time we cut off mesogen–mesogen interactions beyond center-of-mass separations rc = 3.00σ. In addition, we employ a combination of a Verlet and link-cell neighborlist. A mesogen is considered a neighbor of a reference mesogen if their centers of mass are separated by less than rN = 3.50σ.

As we are not interested in simulating any particular material we express all quantities in dimensionless (i.e., “reduced”) units. For example, energy is given in units of εmm, length in units of σ, and temperature in units of εmm/kB. Other derived units are obtained as suitable combinations of these basic ones (see Appendix B.1 of the book by Allen and Tildesley37).

We consider a thermodynamic state characterized by T = 0.95 and P = 1.80 corresponding to a mean number density ρ ≃ 0.92 for which the host phase is sufficiently deep in the nematic phase indicated by Sb ≃ 0.70. For all the simulations and regardless of the spatial arrangement of the colloidal pair we take sz = 24.0 so that the immediate environment of the colloids is not affected directly by the presence of the solid substrates whose sole purpose is to fix [n with combining circumflex]0.

4.2 Bulk phase

We begin our presentation of properties of the bulk phase in the absence of any colloidal particle. To illustrate that under the thermodynamic conditions chosen the host fluid is in the nematic phase we present in Fig. 1 a plot of Sb across the IN phase transition. Because of the relatively small system size the IN phase transition appears to be rounded despite its in principle discontinuous character.
image file: c5sm01860a-f1.tif
Fig. 1 Plot of the nematic order parameter Sb as a function of temperature T (image file: c5sm01860a-u1.tif). Data are shown for N = 1000 mesogens. The IN phase transition occurs at TIN and is obtained from an analysis of the second-order cumulant (see text). The full line represents a spline fitted to the discrete data points and intended to guide the eye. Also shown is SIN = ⅓ at TIN predicted by LdG theory (image file: c5sm01860a-u2.tif) (see text).

Another well-known finite-size effect is visible in the isotropic phase (i.e., for T > TIN) where Sb approaches a small nonzero plateau value of about 0.069. This can be explained by assuming that the liquid crystal consists of molecular-size domains in which mesogens align their longer axes preferentially because of the form of φmm. In the isotropic phase these domains are uncorrelated. However, their number is finite so that by averaging Sb over the ordered domains a residual nonzero value is obtained.

At this stage it is noteworthy that finite system size affects the nematic order parameter only in the isotropic phase and near the IN phase transition [see, for example, Fig. 6(b) of ref. 12 or Fig. 2(b) of ref. 31] whereas Sb is insensitive to system size sufficiently deep in the nematic phase. Hence, under the present thermodynamic conditions (see Section 4.1) a significant system-size effect is not anticipated for the pure host phase (i.e., in the absence of colloidal particles). The presence of the colloids will cause formation of nearly isotropic domains of a certain size determined by the surface curvature of the colloidal particle (i.e., by its size). In these domains, S(r) is small but nonzero. However, size and shape of the domains and the actual value of S(r) reflect the true physics of the system and should not be confused with finite-size effects in the bulk and in the absence of the colloids as discussed before.

Another feature illustrated by the plot in Fig. 1 is that the value predicted by LdG theory

image file: c5sm01860a-t20.tif(26)
agrees remarkably well with SIN ≃ 0.36 obtained from MC using FSS (see below). Moreover, SIN from LdG theory turns out to be universal in that it neither depends on ρ nor T. This is similar to Maier–Saupe theory where a similar universal value is found.38 However, more recent MMF DFT calculations showed that instead SIN exhibits a temperature dependence such that a limiting value of SIN is approached from above if TIN is sufficiently high.12

To estimate TIN (and therefore SIN) we resort to FSS theory and compute g2 from eqn (25) using image file: c5sm01860a-t21.tif for three system sizes corresponding to N = 500, N = 1000, and N = 5000 mesogens. That data obtained for these system sizes are significant and meaningful is concluded from the much more detailed analysis of the IN phase transition in the present model conducted earlier by Greschek and Schoen.31

Results of the present calculations displayed in Fig. 2 indicate that above TIN (i.e., in the isotropic phase), g2 turns out to be independent of system size as one would expect according to the scaling behavior of 〈Sb〉 ∝ N−1/2 and that of 〈Sb2〉 ∝ N−1.26,31 As one approaches TIN from above, g2 first passes through a maximum if N is sufficiently large and then declines sharply with decreasing T where in the nematic phase g2 turns out to be the smaller the larger N is. Most importantly, however, all three curves plotted in Fig. 2 pass through a common intersection demarcating TIN ≃ 1.02 according to the discussion in Section 3.2.3.

image file: c5sm01860a-f2.tif
Fig. 2 Plots of second-order cumulants g2 as functions of temperature T for N = 500 (image file: c5sm01860a-u3.tif), N = 1000 (image file: c5sm01860a-u4.tif), and N = 5000 (image file: c5sm01860a-u5.tif). Inset is an enlargement in the vicinity of the IN phase transition.

Equipped with this result and with expressions for a, B, and C [see eqn (22a) and (22b)], we are now in a position to estimate T* where we find T*/TIN ≃ 0.746. For MBBA, Senbetu and Woo's experimental data allow us to estimate T*/TIN ≃ 0.904 which is of about the same order of magnitude.39 Thus, we conclude that our combined MMF DFT-LdG-FSS approach provides a sufficiently realistic description of the nematic host phase.

With the parameters T*, a, B, and C we plot the LdG free energy density Δf0 from eqn (13) in Fig. 3. As expected, the absolute minimum of Δf0 corresponds to the isotropic phase (Sb = 0) for T > TIN. Exactly at T = TIN, ΔfLdG exhibits a double minimum, one at Sb = 0, the other one at Sb > 0 in the nematic phase. The depth of both minima is the same, that is isotropic and nematic phases coexist. At T slightly below TIN the depth of the minimum at Sb > 0 exceeds that at Sb = 0 so that the nematic phase is thermodynamically stable whereas the isotropic phase is only metastable. Finally, at T = T* the plot of Δf0 exhibits a saddle point at Sb = 0 indicating that the isotropic phase is unstable for all TT*.

image file: c5sm01860a-f3.tif
Fig. 3 Plot of the change in LdG free energy density ΔfLdG as a function of the nematic order parameter S for T > TIN (image file: c5sm01860a-u6.tif), T = TIN (image file: c5sm01860a-u7.tif), T < TIN (image file: c5sm01860a-u8.tif), and T = T* (image file: c5sm01860a-u9.tif).

4.3 Core region

Turning now to our composite system in which colloids are immersed into the nematic host phase, we begin by considering a single colloidal particle first. Plots of the local nematic order parameter S(r) and the director field [n with combining circumflex](r) in Fig. 4(a) reveal the formation of a defect near the colloid's north pole and that this defect is of the Boojum topology as anticipated for a locally planar anchoring of mesogens at the colloid's surface. Upon entering this region, S(r) declines sharply so that the defect has a well-defined boundary.
image file: c5sm01860a-f4.tif
Fig. 4 (a) Director field [n with combining circumflex](r) (dashes) and local nematic order parameter S(r) (see attached color bar) projected onto the xy plane. The grey semicircle marks the upper hemisphere of a colloidal particle with part of a Boojum defect topology centered on its north pole. (b) Three-dimensional sketch of the variation of [n with combining circumflex](r) for a hyperbolic hedgehog defect topology. (c) Variation of S(r) (left ordinate, image file: c5sm01860a-u10.tif), the Frank–Oseen contribution fFO(r) [right ordinate, image file: c5sm01860a-u11.tif, see eqn (14)], and the biaxial order parameter ζ(r) (left ordinate, see text, image file: c5sm01860a-u12.tif) as functions of xr0 and z = 0 cutting through the defect core. Notice that the data plotted have been averaged over a strip of width δy = 1.2 centered on y = 0 to get a reasonably smooth variation of all three quantities shown. The vertical solid line marks the radius Rcore of the circular defect core. Red marks have been added to all three parts of the figure to assist the reader in relating them.

As already mentioned in Section 3.1, special precaution has to be taken to treat the contribution of the defect core to Δf(r). Within the defect core the variation of [n with combining circumflex](r) bears a lot of structural similarity with a hyperbolic hedgehog defect in the bulk [see Fig. 4(b)] represented by [n with combining circumflex](r) = (x, y, −z)T where superscript T denotes the transpose.

Moreover, plots in Fig. 4(c) show that outside the defect core both S(r) and fFO(r) vary rather weakly. Hence, the assumption underlying eqn (14), namely the variation of [n with combining circumflex](r) on a length scale exceeding the molecular one, is well justified. Therefore, felfFO outside the defect core to a good approximation [see eqn (14)].

Inside the defect core, however, [n with combining circumflex](r) varies on a molecular scale such that fFO passes through a relatively sharp maximum and then declines sharply as one approaches the center of the core region [see Fig. 4(c)] so that the assumption underlying fFO(r) is no longer justified. To account for the free-energy contribution of the defect core we therefore resort to a procedure suggested earlier by Lubensky et al.20

These authors derive an analytic expression for the free energy of defect cores considering analytical director fields (such as the one for the hyperbolic hedgehog defect) [see eqn (8) of ref. 20]. Because of the plots in Fig. 4(a) and (b) we take half of this free energy and assign it to a spherical core volume [/]πRcore3 to obtain a free-energy density of fcore = K/Rcore2 where Rcore is the radius of a circular Boojum defect core. We determine the size of the core region by cutting through the center of the defect core of an isolated colloid along the x-axis and take as Rcore that value of x at which S(r) drops below SIN for the first time. Using Rcore ≃ 1.80 [see Fig. 4(c)] we obtain fcore ≈ 0.50kBT which is not an unrealistic value as we conclude by comparison with the much more sophisticated DFT study of defect-core free-energy densities of de las Heras et al.22

To approximate the free-energy density of more complex defect topologies to be discussed in Section 4.4 we assume that fcore is the same everywhere in the core region irrespective of the defect topology. Hence, the free energy of the defect core is obtained from the expression

image file: c5sm01860a-t22.tif(27)
where Vcore = {r|S(r) ≤ SIN|} is the total volume of the defect core. In practice, we determine Vcore by partitioning the entire system into small cubes of side length δs = 0.2 and count the number of cubes in which S(r) ≤ SIN. The total volume of all these small cubes is then equal to Vcore. In a similar fashion we compute Δf(r) for small equally sized volume elements and obtain Δ[scr F, script letter F] through a three-dimensional integration of Δf(r) along the x-, y-, and z-axis using a simple trapezoidal rule.

However, it needs to be stressed that this is truly only a rough approximation to the free energy of defect cores even though it is a standard one in the literature.20,33 Perhaps the most significant assumption behind the expression in eqn (27) is that of insensitivity of fcore to variations in the topology of defects as illustrated by plots in Fig. 4. To improve this situation one could in principle invoke the approach of Lubensky et al.20 and try to find an analytic expression for [n with combining circumflex](r) inside the defect core such that for each and every topology observed a different fcore might obtain analytically. However, if and to what an extent this is possible is currently unknown and would require a study in its own right.

Nevertheless, we feel that the assumption of a assigning the same constant fcore regardless of the specific defect topology is not unrealistically crude. This is concluded from the work of Lubensky et al.20 who show that even for rather disparate [n with combining circumflex](r) the free energy of the defect core is more or less the same such that [scr F, script letter F]coreVcore to a reasonable degree.

4.4 The effective interaction potential

Based upon results presented in Sections 4.2 and 4.3 we now focus on the effective interaction between a pair of colloids immersed in the nematic bulk host phase. We begin in Fig. 5 by considering a pair of colloids in contact with each other. Plots (a)–(c) in the upper panel of Fig. 5 indicate that part of the Boojum defects at isolated colloids interact at sufficiently close proximity of these colloids. For example, for an angle θ = 0° between the intercolloidal center of mass distance r12 and the far-field nematic director [n with combining circumflex]0 a toroidal defect structure exists surrounding the point of contact between the colloids. As θ increases this torus is first “ripped apart” and eventually a handle-like structure forms at θ = 90°.
image file: c5sm01860a-f5.tif
Fig. 5 Upper panel (a)–(c) shows plots of the three-dimensional defect topologies of a pair of colloids (grey spheres) immersed in a nematic host fluid for various angles θ between the center-of-mass distance vector r12 and the far-field nematic director [n with combining circumflex]0 given in the plots. Defect regions are colored in red subject to the condition S(r) ≤ ⅓. Plots in the middle panel (d)–(f) show the corresponding local director field (dashes) and the local nematic order parameter (see attached color bar) projected onto the xy plane where grey circles are similar two-dimensional projections of the colloids. Plots of the biaxiality order parameter are shown in the lower panel (g)–(i). In all cases [n with combining circumflex]0·êx = 1.

Plots (d)–(f) in the middle panel of Fig. 5 are projections of [n with combining circumflex](r) onto the xy plane for the same angles as in parts (a)–(c) of the same figure. These plots indicate that [n with combining circumflex](r) = [n with combining circumflex]0 except in the immediate vicinity of the colloidal pair. As one approaches the colloids, [n with combining circumflex](r) is deformed with respect to [n with combining circumflex]0.

Finally, we plot in parts (g)–(i) of Fig. 5 the local biaxiality order parameter ζ(r). We compute this quantity from our MC simulations via the smallest and middle eigenvalue of Q(r) in eqn (17) which can be expressed as λ = −[S(r) + ζ(r)]/2 and λ0 = −[S(r) − ζ(r)]/2, respectively, for a system with biaxial symmetry because TrQ(r) = 0. From Fig. 5(g)–(i) one notices that biaxiality is relatively small and restricted to the immediate vicinity of the colloids' surfaces. We ascribe this to the nanoscopic size of our colloids (i.e., to the large curvature of their surfaces). A comparison with plots in Fig. 5(d)–(g) illustrates the structural similarity between both sets of plots. More specifically, where S(r) is lowest, ζ(r) is largest.

This also offers another route to treating [scr F, script letter F]core at least in principle because the plots in Fig. 5 suggest that inside the defect core Q(r) remains non-singular. Instead of invoking the approximations introduced in Section 4.3 one could therefore try to use directly Q(r) from eqn (17) obtained in the MC simulations and plug it into the right-hand side of the expression on the first line of eqn (14). The required differentiation of Q(r) has, of course, to be performed numerically but would allow one to compute fel(r) also inside the defect core.

Unfortunately, in practice we found that our data are way too noisy to follow this route and obtain reliable values for fel(r) inside the defect core. A reliable numerical differentiation of Q(r) would need a much finer discretization of the grid on which this quantity is stored because of its rather strong spatial variation inside the defect core. Notice, that this does not contradict the smoothness of plots in Fig. 4(c) as the data shown there have been averaged over a fairly wide strip inside the defect core.

The magnitude of the spatial variation of Q(r) is reflected in part by the plot of fFO(r) in Fig. 4(c) sufficiently deep inside the defect core. In this region, fFO(r) exceeds its typical values in the elastic regime just outside this region by about two to three orders of magnitude. As is clear from eqn (14) the large value of fFO(r) can immediately be traced back to a very strong spatial variation of [n with combining circumflex](r) inside the defect core which is linked to a similarly strong spatial variation of Q(r) because of eqn (9).

Consequently, to be able to differentiate this quantity numerically and reliably would require a better resolution of Q(r) on a much finer grid. This obviously would entail much longer MC simulations to obtain reasonably good statistics which, unfortunately, is beyond reach given the number of simulations and the typical system sizes necessary to map out the effective interaction potential between the colloidal pair with sufficient resolution. Because of these constraints and because of the discussion in Section 4.3 we conclude that our present treatment of the free-energy contribution of the defect core offers the best possible approximation.

In addition to a distortion of [n with combining circumflex](r) plots in Fig. 5(d)–(f) reveal regions in which S(r) ≪ Sb. From the parallel three-dimensional plots in Fig. 5(a)–(c) one sees that the volume of these regions changes with θ. Hence, it seems intuitive to introduce an associated change in effective free energy

Δ[scr F, script letter F]eff = [scr F, script letter F]el + [scr F, script letter F]core + Δ[scr F, script letter F]LdG − 2Δ[scr F, script letter F]B(28)
where Δ[scr F, script letter F]B is the change in free energy associated with a single Boojum defect. In computing Δ[scr F, script letter F]B we assume that it also consists of an elastic Frank, a defect-core, and a LdG free energy contribution. In practice, it turns out that the contribution of each of the first three terms on the right-hand side of eqn (28) for the interacting colloidal pair relative to the corresponding contribution to Δ[scr F, script letter F]B is roughly of the same order of magnitude in the range of kBT in agreement with plots presented in Fig. 6 and in Fig. 7.

image file: c5sm01860a-f6.tif
Fig. 6 The effective free energy Δ[scr F, script letter F]eff in units of the change in free energy associated with an isolated Boojum colloid Δ[scr F, script letter F]B. Vertical dashed lines demarcate minima in the curves plotted; the limiting value θ ≈ 49° is also indicated (see text); r12 = 7.20 (image file: c5sm01860a-u13.tif), r12 = 8.00 (image file: c5sm01860a-u14.tif).

image file: c5sm01860a-f7.tif
Fig. 7 Contour plot of Δ[scr F, script letter F]eff[scr F, script letter F]B (see attached color bar) as a function of rT12 = (x12, y12, 0)T. Curves plotted in Fig. 6 correspond to r12 = 7.20 (image file: c5sm01860a-u15.tif) and to r12 = 8.00 (image file: c5sm01860a-u16.tif), respectively.

Data for Δ[scr F, script letter F]eff plotted in Fig. 6 exhibit a number of interesting features. First, one notices that depending on θ, the sum [scr F, script letter F]el + [scr F, script letter F]core + Δ[scr F, script letter F]LdG may be larger or smaller than 2Δ[scr F, script letter F]B and therefore Δ[scr F, script letter F]eff may be viewed as a repulsive or attractive effective interaction potential acting between a colloidal pair and mediated by the nematic host phase. Second, the minimum in the plot of Δ[scr F, script letter F]eff shifts to larger angles θmin as the intercolloidal center-of-mass distance increases. This is fully in line with experimental observations made by Smalyukh et al.11 and theoretical observations made by Tasinkevych et al.8 Third, these latter authors could also demonstrate that their experimental value of θmin increases monotonically towards a limiting value. This limiting value is found by realizing that for a single Boojum defect (corresponding to a sufficiently separated pair of colloids) the spatial variation of the nematic director field bears close resemblance to the spatial variation of the electrostatic field associated with interacting quadrupoles. For the latter the orientation dependence of the electrostatic energy can be cast as

U ∝ 9 − 90[thin space (1/6-em)]cos2[thin space (1/6-em)]θ+ 105[thin space (1/6-em)]cos4[thin space (1/6-em)]θ(29)
from which θmin ≈ 49° follows without further ado.

One also notices from the plots in Fig. 6 an upward shift of the curves that causes larger spatial regions in which the effective potential Δ[scr F, script letter F]eff is repulsive. In particular, such a repulsive region exists for angles θ ≳ 70° and r12 = 8.00. This is fully in line with trajectories of colloids measured by Smalyukh et al.11 by video microscopy and presented in their Fig. 2.

From curves such as the ones presented in Fig. 6 we are now in a position to present in Fig. 7 a more refined contour plot of Δ[scr F, script letter F]eff illustrating in a broader way the structural complexity of the effective-potential landscape.

Using the data plotted in Fig. 7 we can now also investigate structures that several colloids would form in a nematic host phase. For simplicity, and because the experiments used a quasi two-dimensional setup11 we consider a two-dimensional, coarse-grained system treating the nematic host implicitly via Δ[scr F, script letter F]eff. To account for the evolution of the colloids in configuration space we employ a conventional canonical-ensemble Metropolis MC scheme.

To that end we store Δ[scr F, script letter F]eff on a two-dimensional regular grid with a spacing of 0.2 between neighboring nodes. As a position of a colloid will normally not coincide with a grid point, we interpolate Δ[scr F, script letter F]eff between the four nearest-neighbor grid points in a bilinear fashion to get Δ[scr F, script letter F]eff at the actual center-of-mass position of a colloidal disk. A displacement of a disk is then accepted with a probability min[1, exp(−Δ[scr F, script letter F]eff[scr F, script letter F]B)].

For two packing fractions ϕ = Ncollπr02/sxsy of Ncoll colloidal disks, characteristic “snapshots” from the simulations at thermodynamic equilibrium are shown in Fig. 8. Here sx and sy are linear dimensions of the simulation cell. As one can see colloidal disks tend to form linear chains at the lower packing fraction where the symmetry axis of each chain forms an angle of θ ≈ 30° with [n with combining circumflex]0 [see Fig. 8(a)] whereas at higher packing fraction a more extended two-dimensional network of colloidal disks exists at thermodynamic equilibrium [see Fig. 8(b)]. Both types of structures bear a remarkable resemblance to those displayed in Fig. 1(b) and (e) of the work by Smalyukh et al.11

image file: c5sm01860a-f8.tif
Fig. 8 “Snapshots” from canonical-ensemble MC simulations of colloidal disks immersed in a nematic host phase taken into account implicitly via the effective interaction potential Δ[scr F, script letter F]eff shown in Fig. 7. (a) ϕ = 0.065, (b) ϕ = 0.234 (sx = sy = 50 and [n with combining circumflex]0·êx = 1). The direction of the far-field director [n with combining circumflex]0 is indicated in the figure.

5 Discussion and conclusions

By means of a novel hybrid approach we investigate the self-assembly of spherical colloidal particles with chemically homogeneous surfaces at a coarse-grained level of description. Our hybrid approach combines molecular-scale methods and theories such as MC computer simulations, classical DFT, and elements of FSS theory with macroscopic theories such as LdG and the Frank–Oseen treatment of the free-energy density associated with elastic deformations of the director field. The goal of this approach is to integrate out less relevant degrees of freedom in a controlled fashion while maintaining the correct physics of a composite system such as the present one. The philosophy of our approach is quite general and may perhaps be applied to other composite systems as well. An example in this respect could be that of Janus colloids immersed into a nematic host phase.

In our case, the colloids are immersed into a nematic liquid-crystal host phase. Self-assembly is driven by effective interactions mediated by the host. The effective interactions correspond to a change in free energy caused by defects (i.e., regions of lower nematic order) in the host phase. These defects cause a local decline of nematic order and distortion of the director field. Both features arise because of the mismatch between the local alignment of mesogens at the curved surface of the colloids and the nematic far-field director [n with combining circumflex]0. They depend on the arrangement of the colloids in space.

We are treating both the local lowering of nematic order and the distortion of the director field at mean-field level assuming that the former can be adequately accounted for by LdG theory and the latter within the standard expression for the Frank elastic free energy. The mean-field character is reflected by the fact that we compute LdG expansion parameters a, B, and C from mean-field DFT and that the elastic contribution to the change in free energy does not account for fluctuations of [n with combining circumflex](r). However, the Frank free energy is computed only outside the defect core where [n with combining circumflex](r) changes on a length scale large compared to a molecular one.

Applying LdG theory within the framework of molecular simulations poses a fundamental problem because under most conditions statistical accuracy is insufficient to compute the LdG expansion coefficients B and C.26 However, employing classical DFT allows one to express analytically the LDG coefficients in terms of thermodynamic state parameters such as density and temperature.

A key ingredient in LdG theory is the temperature T* at which the isotropic phase becomes thermodynamically unstable. It is related to the temperature TIN at which one has IN phase coexistence. To improve our DFT estimate of T* for the LdG treatment we locate TIN through an analysis of second-order cumulants of the nematic order parameter within the framework of FSS theory.

An alternative route to the LdG coefficients has recently been proposed by Gupta and Ilg.40 Starting from a state point in the isotropic phase for which they can compute the free energy semi-analytically, Gupta and Ilg apply an external ordering field to drive the system into an ordered nematic phase. In the ordered state the free energy can then be calculated by thermodynamic integration. The intrinsic free-energy contribution can be related to functions that can, on the one hand, be determined numerically through reliable fit functions and that one can, on the other hand relate to the LdG expansion coefficients analytically.

A crucial ingredient in the approach of Gupta and Ilg is the ordering field which for their Gay–Berne model of a liquid crystal has a magnitude of the order of one. When applying the technique of Gupta and Ilg to our model system we found that the magnitude of the ordering field had to be about five orders of magnitude larger to drive a system from the isotropic to the nematic phase. We believe that this huge difference in the external field is caused by the almost spherical shape of mesogens in our model. Because of the disparate magnitude of the ordering fields it turned out that for the present model the approach of Gupta and Ilg could not be applied reliably.

However, our combined MC-DFT-FSS approach allows us to compute the effective interaction potential reliably. A comparison with experimental data reveals that

(1) the effective potential has a minimum at an angle θ ≈ 30° between the intercolloidal distance vector r12 and the far-field nematic director [n with combining circumflex]0 if the colloids are sufficiently close to each other.

(2) the position of the minimum shifts monotonically to larger θ if r12 increases.

(3) depending on θ and r12 the effective potential may be repulsive or attractive.

These features turn out to be in semi-quantitative agreement with experimental observations.9,11

Encouraged by the apparent consistency of these observations with experimental data we perform coarse-grained, canonical-ensemble MC simulations of several colloidal disks immersed in a nematic host phase that is now treated implicitly through the effective interaction potential. In that regard our approach has a multiscale character. The structures observed at different packing fractions agree again qualitatively with experimental data11 despite the much larger colloids used experimentally.

However, it needs to be stressed that the simulations of several colloids is based upon the assumption of pairwise additivity of the effective interactions. A priori there is no guarantee that this assumption is valid. However, the excellent qualitative agreement between structures observed within our approach with those seen experimentalls seems to justfy the assumption of pairwise additive effective interactions a posteriori.

In summary, we have shown that the self-assembly of colloidal particles in a nematic liquid crystal is driven by the occurrence of colloid-induced defects associated with a local elastic distortion of the director field. This also implies that in an isotropic phase no such self-assembly would occur as the effective interaction potential would vanish identically.


We are grateful for financial support from the Deutsche Forschungsgemeinschaft via the International Graduate Research Training Group 1524.


  1. P. G. de Gennes and J. Prost, The physics of liquid crystals, Clarendon, Oxford, 2nd edn, 1995 Search PubMed.
  2. V. S. R. Jampani, M. Škarabot, M. Ravnik, S. Čopar, S. Žumer and I. Muševič, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 031703 CrossRef CAS PubMed.
  3. K. Izaki and Y. Kimura, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 062507 CrossRef PubMed.
  4. U. Ognysta, A. Nych, V. Nazarenko, M. Škarabot and I. Muševič, Langmuir, 2009, 25, 12092–12100 CrossRef CAS PubMed.
  5. H. Qi and T. Hegmann, J. Mater. Chem., 2006, 16, 4197–4205 RSC.
  6. M. Humar, M. Ravnik, S. Pajk and I. Muševič, Nat. Photonics, 2009, 3, 595–600 CrossRef CAS.
  7. I. Muševič, M. Škarabot and M. Humar, J. Phys.: Condens. Matter, 2011, 23, 284112 CrossRef PubMed.
  8. M. Tasinkevych, N. M. Silvestre and M. M. Telo da Gama, New J. Phys., 2012, 14, 073030 CrossRef.
  9. P. Poulin and D. Weitz, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 626–637 CrossRef CAS.
  10. P. Poulin, H. Stark, T. C. Lubensky and D. A. Weitz, Science, 1997, 275, 1770–1773 CrossRef CAS PubMed.
  11. I. I. Smalyukh, O. D. Lavrentovich, A. N. Kuzmin, A. V. Kachynski and P. N. Prasad, Phys. Rev. Lett., 2005, 95, 157801 CrossRef CAS PubMed.
  12. S. Giura and M. Schoen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 022507 CrossRef PubMed.
  13. W. Maier and A. Saupe, Z. Naturforsch., A: Phys. Sci., 1958, 13, 564–566 Search PubMed.
  14. W. Maier and A. Saupe, Z. Naturforsch., A: Phys. Sci., 1959, 14, 882–889 Search PubMed.
  15. B. Jérôme, Rep. Prog. Phys., 1991, 54, 391–451 CrossRef.
  16. M. Škarabot, M. Ravnik, S. Žumer, U. Tkalec, I. Poberaj, D. Babič, N. Osterman and I. Muševič, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 031705 CrossRef PubMed.
  17. A. A. Sonin, The Surface Physics of Liquid Crystals, Gordon and Breach, Amsterdam, 1995 Search PubMed.
  18. M. Melle, S. Schlotthauer, M. G. Mazza, S. H. L. Klapp and M. Schoen, J. Chem. Phys., 2012, 136, 194703 CrossRef PubMed.
  19. R. W. Ruhwandl and E. M. Terentjev, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 5561–5566 CrossRef CAS.
  20. T. C. Lubensky, D. Pettey, N. Currier and H. Stark, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 610–625 CrossRef CAS.
  21. S. Grollau, N. L. Abbott and J. J. de Pablo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 011702 CrossRef CAS PubMed.
  22. D. de las Heras, L. Mederos and E. Velasco, Liq. Cryst., 2010, 37, 45–56 CrossRef CAS.
  23. I. Bajc, F. Hecht and S. Žumer, 2015, arXiv:1505.07046v1.
  24. J. Dontabhaktuni, M. Ravnik and S. Žumer, Soft Matter, 2012, 8, 1657–1663 RSC.
  25. I. Muševič, U. Škarabot, U. Tkalec, M. Ravnik and S. Žumer, Science, 2006, 313, 954–958 CrossRef PubMed.
  26. R. Eppenga and D. Frenkel, Mol. Phys., 1984, 52, 1303–1334 CrossRef CAS.
  27. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical recipes, Cambridge University Press, Cambridge, 3rd edn, 2007 Search PubMed.
  28. M. P. Allen and D. Frenkel, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 37, 1813R–1816R CrossRef.
  29. M. P. Allen and D. Frenkel, Phys. Rev. A: At., Mol., Opt. Phys., 1990, 42, 3641 CrossRef.
  30. T. Stieger, M. Schoen and M. G. Mazza, J. Chem. Phys., 2014, 140, 054905 CrossRef PubMed.
  31. M. Greschek and M. Schoen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 011704 CrossRef PubMed.
  32. S. Giura, B. G. Márkus, S. H. L. Klapp and M. Schoen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 012313 CrossRef PubMed.
  33. M. Kleman and O. D. Lavrentovich, Soft Matter Physics: An introduction, Springer-Verlag, New York, 2003 Search PubMed.
  34. H. Weber, W. Paul and K. Binder, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 2168–2174 CrossRef CAS.
  35. K. Vollmayer, J. D. Reger, M. Scheucher and K. Binder, Z. Phys. B: Condens. Matter, 1993, 91, 113–125 CrossRef.
  36. D. Frenkel and B. Smit, Understanding molecular simulation: From algorithms to applications, Academic Press, San Diego, 2nd edn, 2002 Search PubMed.
  37. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, Oxford, 2002 Search PubMed.
  38. W. Maier and A. Saupe, Z. Naturforsch., A: Phys. Sci., 1960, 15, 287–292 Search PubMed.
  39. L. Senbetu and C.-W. Woo, Mol. Cryst. Liq. Cryst., 1982, 84, 101–124 CrossRef CAS.
  40. B. Gupta and P. Ilg, Polymers, 2013, 5, 328–343 CrossRef.

This journal is © The Royal Society of Chemistry 2016