Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Sergej
Püschel-Schlotthauer
^{a},
Tillmann
Stieger
^{a},
Michael
Melle
^{a},
Marco G.
Mazza
^{b} and
Martin
Schoen
^{ac}
^{a}Stranski-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
^{b}Max-Planck-Institut für Dynamik und Selbstorganisation, Am Faßberg 17, 37077 Göttingen, Germany
^{c}Department 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.

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 Weitz^{9} found experimentally that in a nematic phase the colloidal center-to-center distance vector r_{12} forms a “magic” angle of θ ≈ 30° with the global director _{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 Weitz^{9} 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 r_{12} = |r_{12}|; this minimum shifts to larger θ as r_{12} 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.

φ_{mm}(r_{ij}, ω_{i}, ω_{j}) = φ_{iso}(r_{ij}) + φ_{anis}(r_{ij}, ω_{i}, ω_{j}) | (1) |

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

(2) |

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}(r_{ij}, ω_{i}, ω_{j}) = φ_{att}(r_{ij})Ψ(_{ij}, ω_{i}, ω_{j}) | (3) |

Ψ(_{ij}, ω_{i}, ω_{j}) = 5ε_{1}P_{2}[û(ω_{i})·û(ω_{j})] + 5ε_{2}{P_{2}[û(ω_{i})·_{ij}] + P_{2}[û(ω_{j})·_{ij}]} | (4) |

To predict and control _{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 _{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 s_{z} along the z-axis. Thus, the interaction between a mesogen and the substrates can be cast as

(5) |

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 _{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) |

(7) |

In eqn (7), η is the inverse Debye screening length and parameters a_{1} and a_{2} have been introduced to guarantee that the minimum of φ_{mc} remains at a distance r_{0} + σ 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

g_{c}(^{c}_{ij}, ω_{i}) = [1 − |^{c}_{ij}·û(ω_{i})|]^{2} | (8) |

(9) |

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

Δf(r) = Δf_{LdG}(r) + f_{el}(r) + f_{core} | (10) |

(11) |

(12) |

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

Δf_{0} = AS_{b}^{2} + BS_{b}^{3} + CS_{b}^{4} | (13) |

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

(14) |

(15) |

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 f_{core} 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 functional^{23}

(16) |

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) = S_{b} and (r) = _{0}, Δf_{LdG} = f_{el} = 0.

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 (r) for a thermodynamic equilibrium situation via suitably defined ensemble averages. One may then feed S(r) and (r) into expressions such as eqn (12) and (14) to eventually obtain the absolute minimum of Δ [after an integration of Δf(r) over volume, see also Section 4.3].

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 (r).

Thus, to eventually compute Δ we begin by introducing the local alignment tensor

(17) |

However, to compute Δf_{LdG}(r) and f_{el}(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 (k) of Q(r). In the limit of |k| → 0 simple linear relationships exist from which the Frank constants K_{1}, K_{2}, and K_{3} associated with bend, twist, and splay deformations of (r) can be estimated reliably. Stieger et al.^{30} have recently applied the method of Allen and Frenkel^{28,29} to show that for the present model of the host phase K_{1} ≃ K_{2} ≃ K_{3} = K so that the one-constant form^{1} of f_{el} 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(S_{b}) obtained for a bulk nematic phase without colloidal particles. However, as discussed in detail by Eppenga and Frenkel^{26} and later by Greschek and Schoen^{31} it is next to impossible to determine B and C with sufficient statistical accuracy from P(S_{b}).

where β = 1/k_{B}T (k_{B} is Boltzmann's constant), x = cosθ, and θ is the azimuthal angle if we take the z-axis of our coordinate system to coincide with _{0}. Members of the set {S_{l}} defined on the interval [0, 1] are order parameters and {u_{l}} are parameters that account for the contribution of anisotropic mesogen–mesogen interactions to the free energy, respectively.^{32}

which implies that in the isotropic phase (x) = ½. Again, because of the uniaxiality of the nematic phase we expand (x) in terms of Legendre polynomials {P_{l}} via

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 (ξ^{4}) allows us to recast eqn (18) as

if we limit ourselves to the leading term l = 2 in the expression for ξ and use S_{2} = S_{b}. In eqn (21), a = 2ρk_{B}/5 and T* = −5ρu_{2}/2k_{B} where the latter is the temperature at which the isotropic phase becomes thermodynamically unstable. Assuming Δf_{n} = Δf_{0} one easily obtains

by comparing terms of equal power in S_{b} 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}

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

(18) |

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

(19) |

(20) |

(21) |

(22a) |

(22b) |

3.2.3 Elements of finite-size scaling theory.
Unfortunately, the expression for T* given in the preceding section depends on u_{2} 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, u_{2} = − 32πε_{1}ε_{mm}/15. A more elaborate, temperature dependent form for u_{2} 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.

Second, with an improved estimate for the temperature T_{IN} 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.

of the order-parameter distribution the key quantity in FSS are cumulants of its various moments. Here, the second-order cumulant

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 intersection^{31} which scales as ^{−d} where is the linear extent of a system of dimension d.^{34}

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 that^{33}

(23) |

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

(24) |

(25) |

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 ^{−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.

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 × 10^{5} MC cycles for equilibration followed by another 1.0 × 10^{6} cycles during which relevant ensemble averages are taken. To save computer time we cut off mesogen–mesogen interactions beyond center-of-mass separations r_{c} = 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 r_{N} = 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}/k_{B}. Other derived units are obtained as suitable combinations of these basic ones (see Appendix B.1 of the book by Allen and Tildesley^{37}).

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 S_{b} ≃ 0.70. For all the simulations and regardless of the spatial arrangement of the colloidal pair we take s_{z} = 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 _{0}.

Another well-known finite-size effect is visible in the isotropic phase (i.e., for T > T_{IN}) where S_{b} 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 S_{b} 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 S_{b} 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

(26) |

To estimate T_{IN} (and therefore S_{IN}) we resort to FSS theory and compute g_{2} from eqn (25) using 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 T_{IN} (i.e., in the isotropic phase), g_{2} turns out to be independent of system size as one would expect according to the scaling behavior of 〈S_{b}〉 ∝ N^{−1/2} and that of 〈S_{b}^{2}〉 ∝ N^{−1}.^{26,31} As one approaches T_{IN} from above, g_{2} first passes through a maximum if N is sufficiently large and then declines sharply with decreasing T where in the nematic phase g_{2} 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 T_{IN} ≃ 1.02 according to the discussion in Section 3.2.3.

Fig. 2 Plots of second-order cumulants g_{2} as functions of temperature T for N = 500 (), N = 1000 (), and N = 5000 (). 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*/T_{IN} ≃ 0.746. For MBBA, Senbetu and Woo's experimental data allow us to estimate T*/T_{IN} ≃ 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 Δf_{0} from eqn (13) in Fig. 3. As expected, the absolute minimum of Δf_{0} corresponds to the isotropic phase (S_{b} = 0) for T > T_{IN}. Exactly at T = T_{IN}, Δf_{LdG} exhibits a double minimum, one at S_{b} = 0, the other one at S_{b} > 0 in the nematic phase. The depth of both minima is the same, that is isotropic and nematic phases coexist. At T slightly below T_{IN} the depth of the minimum at S_{b} > 0 exceeds that at S_{b} = 0 so that the nematic phase is thermodynamically stable whereas the isotropic phase is only metastable. Finally, at T = T* the plot of Δf_{0} exhibits a saddle point at S_{b} = 0 indicating that the isotropic phase is unstable for all T ≤ T*.

Fig. 3 Plot of the change in LdG free energy density Δf_{LdG} as a function of the nematic order parameter S for T > T_{IN} (), T = T_{IN} (), T < T_{IN} (), and T = T* (). |

Fig. 4 (a) Director field (r) (dashes) and local nematic order parameter S(r) (see attached color bar) projected onto the x–y 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 (r) for a hyperbolic hedgehog defect topology. (c) Variation of S(r) (left ordinate, ), the Frank–Oseen contribution f_{FO}(r) [right ordinate, , see eqn (14)], and the biaxial order parameter ζ(r) (left ordinate, see text, ) as functions of x − r_{0} 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 R_{core} 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 (r) bears a lot of structural similarity with a hyperbolic hedgehog defect in the bulk [see Fig. 4(b)] represented by (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 f_{FO}(r) vary rather weakly. Hence, the assumption underlying eqn (14), namely the variation of (r) on a length scale exceeding the molecular one, is well justified. Therefore, f_{el} ≈ f_{FO} outside the defect core to a good approximation [see eqn (14)].

Inside the defect core, however, (r) varies on a molecular scale such that f_{FO} 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 f_{FO}(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 πR_{core}^{3} to obtain a free-energy density of f_{core} = K/R_{core}^{2} where R_{core} 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 R_{core} that value of x at which S(r) drops below S_{IN} for the first time. Using R_{core} ≃ 1.80 [see Fig. 4(c)] we obtain f_{core} ≈ 0.50k_{B}T 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 f_{core} 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

(27) |

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 f_{core} 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 (r) inside the defect core such that for each and every topology observed a different f_{core} 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 f_{core} 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 (r) the free energy of the defect core is more or less the same such that _{core} ∝ V_{core} to a reasonable degree.

Plots (d)–(f) in the middle panel of Fig. 5 are projections of (r) onto the x–y plane for the same angles as in parts (a)–(c) of the same figure. These plots indicate that (r) = _{0} except in the immediate vicinity of the colloidal pair. As one approaches the colloids, (r) is deformed with respect to _{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 _{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 f_{el}(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 f_{el}(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 f_{FO}(r) in Fig. 4(c) sufficiently deep inside the defect core. In this region, f_{FO}(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 f_{FO}(r) can immediately be traced back to a very strong spatial variation of (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 (r) plots in Fig. 5(d)–(f) reveal regions in which S(r) ≪ S_{b}. 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

Δ_{eff} = _{el} + _{core} + Δ_{LdG} − 2Δ_{B} | (28) |

Fig. 7 Contour plot of Δ_{eff}/Δ_{B} (see attached color bar) as a function of r^{T}_{12} = (x_{12}, y_{12}, 0)^{T}. Curves plotted in Fig. 6 correspond to r_{12} = 7.20 () and to r_{12} = 8.00 (), respectively. |

Data for Δ_{eff} plotted in Fig. 6 exhibit a number of interesting features. First, one notices that depending on θ, the sum _{el} + _{core} + Δ_{LdG} may be larger or smaller than 2Δ_{B} and therefore Δ_{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 Δ_{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 − 90cos^{2}θ+ 105cos^{4}θ | (29) |

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 Δ_{eff} is repulsive. In particular, such a repulsive region exists for angles θ ≳ 70° and r_{12} = 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 Δ_{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 setup^{11} we consider a two-dimensional, coarse-grained system treating the nematic host implicitly via Δ_{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 Δ_{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 Δ_{eff} between the four nearest-neighbor grid points in a bilinear fashion to get Δ_{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(−Δ_{eff}/Δ_{B})].

For two packing fractions ϕ = N_{coll}πr_{0}^{2}/s_{x}s_{y} of N_{coll} colloidal disks, characteristic “snapshots” from the simulations at thermodynamic equilibrium are shown in Fig. 8. Here s_{x} and s_{y} 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 _{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}

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 Δ_{eff} shown in Fig. 7. (a) ϕ = 0.065, (b) ϕ = 0.234 (s_{x} = s_{y} = 50 and _{0}·ê_{x} = 1). The direction of the far-field director _{0} is indicated in the figure. |

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 _{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 (r). However, the Frank free energy is computed only outside the defect core where (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 T_{IN} at which one has IN phase coexistence. To improve our DFT estimate of T* for the LdG treatment we locate T_{IN} 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 r_{12} and the far-field nematic director _{0} if the colloids are sufficiently close to each other.

(2) the position of the minimum shifts monotonically to larger θ if r_{12} increases.

(3) depending on θ and r_{12} 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 data^{11} 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.

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

This journal is © The Royal Society of Chemistry 2016 |