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

Abstract

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.¹ 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.²

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.³ 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.⁸ 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⁹ found experimentally that in a nematic phase the colloidal center-to-center distance vector r₁₂ forms a “magic” angle of θ ≈ 30° with the global director n̂₀ if the liquid crystal molecules (mesogens) favor a local planar anchoring at the colloid's surface.¹⁰

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.⁹ 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”.⁹

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

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 θ.¹¹ 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̂₀. 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(r_ij, ω_i, ω_j) = φ_iso(r_ij) + φ_anis(r_ij, ω_i, ω_j) (1)

where r_ij = r_i − r_j is the distance vector connecting the centers of mass of a mesogenic pair located at points r_i and r_j, respectively, and r_ij = |r_ij|. 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

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.¹² From a lengthy but relatively straightforward derivation these authors show that

φ_anis(r_ij, ω_i, ω_j) = φ_att(r_ij)Ψ(r̂_ij, ω_i, ω_j) (3)

In eqn (3) the anisotropy function is given by

Ψ(r̂_ij, ω_i, ω_j) = 5ε₁P₂[û(ω_i)·û(ω_j)] + 5ε₂{P₂[û(ω_i)·r̂_ij] + P₂[û(ω_j)·r̂_ij]} (4)

where û(ω_i,j) and r̂_ij = r_ij/r_ij 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, P₂(x) = ½(3x² − 1) is the second Legendre polynomial and the (dimensionless) anisotropy parameters 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 model^13,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.¹² Unfortunately, when a nematic phase is forming in the bulk it is impossible to determine beforehand the direction of n̂₀. In fact, there is an infinite number of easy axes¹⁵ with which n̂₀ may align in the bulk nematic phase.

To predict and control n̂₀ 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.¹⁶ Because orientational order in a nematic liquid crystal is a long-range phenomenon, substrate-induced alignment of mesogens allows one to control n̂₀ 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.¹⁷

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

where Δz_i = z_i ± s_z/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̂₀ 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]² (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̂₀ to the x-axis on average where |n̂₀·ê_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 r₀ = 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 us¹⁸ we take the mesogen–colloid interaction potential to be given bywhere r^c_ij = |r_i − r^c_j| is the distance between the center of mass of mesogen i at r_i and that of colloid j at r^c_j. Hence, r^c_ij − r₀ 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 a₁ and a₂ have been introduced to guarantee that the minimum of φ_mc remains at a distance r₀ + σ from the colloid's center of mass and to maintain the depth of the attractive well at ε_mc irrespective of η.¹⁸ 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(r̂^c_ij, ω_i) = [1 − |r̂^c_ij·û(ω_i)|]² (8)

where r̂^c_ij = r^c_ij/r^c_ij is the local surface normal of the colloid. Hence, the anchoring function introduced in eqn (8) serves to align mesogens in a degenerate,¹⁵ 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̂(r) to deviate from n̂₀ in certain regions near the colloid's surface. Associated with this distortion of n̂(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̂(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 aswhere S(r) is the nematic order parameter, n_α(r) is the α-component of n̂(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) = Δf_LdG(r) + f_el(r) + f_core (10)

where the Landau–de Gennes contribution is given byusing 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_γα = ¾S³ This allows us to rewrite eqn (11) as 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)

Δf₀ = AS_b² + BS_b³ + CS_b⁴ (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), S_b 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

where f_FO(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 To apply eqn (14) some caution is advisable. This is because the expression for f_FO(r) in eqn (14) is derived under the assumption that spatial variations of n̂(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 f_FO(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) ≤ S_IN where S_IN 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.²² 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²³

However, this procedure has a twofold drawback. First, the minimization bears the danger that its solutions S(r) and n̂(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̂(r) in thermodynamic equilibrium is rather complex.²⁴ 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.²⁵

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 n̂(r) = n̂₀, Δf_LdG = f_el = 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̂(r) for a thermodynamic equilibrium situation via suitably defined ensemble averages. One may then feed S(r) and n̂(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

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

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 Q̃(k) of Q(r). In the limit of |k| → 0 simple linear relationships exist from which the Frank constants K₁, K₂, and K₃ associated with bend, twist, and splay deformations of n̂(r) can be estimated reliably. Stieger et al.³⁰ have recently applied the method of Allen and Frenkel^28,29 to show that for the present model of the host phase K₁ ≃ K₂ ≃ K₃ = K so that the one-constant form¹ of f_el in eqn (14) is an excellent approximation. Under the thermodynamic conditions used here (see Section 4.1), K ≃ 1.66ε_mmσ⁻¹ 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²⁶ and later by Greschek and Schoen³¹ it is next to impossible to determine B and C with sufficient statistical accuracy from P(S_b).

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

where β = 1/k_BT (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 n̂₀. 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.³²

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

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 {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 [scr O, script letter O] (ξ⁴) allows us to recast eqn (18) asif we limit ourselves to the leading term l = 2 in the expression for ξ and use S₂ = S_b. In eqn (21), a = 2ρk_B/5 and T* = −5ρu₂/2k_B where the latter is the temperature at which the isotropic phase becomes thermodynamically unstable. Assuming Δf_n = Δf₀ one easily obtainsby 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.¹

3.2.3 Elements of finite-size scaling theory. Unfortunately, the expression for T* given in the preceding section depends on u₂ 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₂ = − 32πε₁ε_mm/15. A more elaborate, temperature dependent form for u₂ 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.³¹ This failure can be linked to the complete neglect of pair correlations at SMF level and their overestimation in liquidlike phases at MMF level.¹²

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³³

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.

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

of the order-parameter distribution the key quantity in FSS are cumulants of its various moments. Here, the second-order cumulantis particularly useful.³¹ 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³¹ which scales as [script L]^−d where [script L] is the linear extent of a system of dimension d.³⁴

Moreover, it has been demonstrated by Vollmayer et al.³⁵ 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.¹⁸ 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 s_x/s_y, the distance s_z between the solid substrates, and the transverse component P_∥ = ½(P_xx − P_yy) 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.³⁶

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⁵ MC cycles for equilibration followed by another 1.0 × 10⁶ 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³⁷).

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 n̂₀.

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 S_b 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.

Fig. 1 Plot of the nematic order parameter S_b as a function of temperature T (). Data are shown for N = 1000 mesogens. The IN phase transition occurs at T_IN 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 S_IN = ⅓ at T_IN predicted by LdG theory () (see text).

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

agrees remarkably well with S_IN ≃ 0.36 obtained from MC using FSS (see below). Moreover, S_IN 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.³⁸ However, more recent MMF DFT calculations showed that instead S_IN exhibits a temperature dependence such that a limiting value of S_IN is approached from above if T_IN is sufficiently high.¹²

To estimate T_IN (and therefore S_IN) we resort to FSS theory and compute g₂ 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.³¹

Results of the present calculations displayed in Fig. 2 indicate that above T_IN (i.e., in the isotropic phase), g₂ 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²〉 ∝ N⁻¹.^26,31 As one approaches T_IN from above, g₂ first passes through a maximum if N is sufficiently large and then declines sharply with decreasing T where in the nematic phase g₂ 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₂ 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.³⁹ 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₀ from eqn (13) in Fig. 3. As expected, the absolute minimum of Δf₀ 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₀ 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* ().

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̂(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.

Fig. 4 (a) Director field n̂(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 n̂(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₀ 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 n̂(r) bears a lot of structural similarity with a hyperbolic hedgehog defect in the bulk [see Fig. 4(b)] represented by n̂(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 n̂(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, n̂(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.²⁰

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³ to obtain a free-energy density of f_core = K/R_core² 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_BT 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.²²

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

where V_core = {r|S(r) ≤ S_IN|} is the total volume of the defect core. In practice, we determine V_core by partitioning the entire system into small cubes of side length δs = 0.2 and count the number of cubes in which S(r) ≤ S_IN. The total volume of all these small cubes is then equal to V_core. 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 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.²⁰ and try to find an analytic expression for n̂(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.²⁰ who show that even for rather disparate n̂(r) the free energy of the defect core is more or less the same such that [scr F, script letter F]_core ∝ V_core 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 r₁₂ and the far-field nematic director n̂₀ 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°.

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 r₁₂ and the far-field nematic director n̂₀ 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 x–y 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̂₀·ê_x = 1.

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

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 λ₀ = −[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 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 n̂(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̂(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

Δ[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 k_BT in agreement with plots presented in Fig. 6 and in Fig. 7.

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); r₁₂ = 7.20 (), r₁₂ = 8.00 ().

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 r^T₁₂ = (x₁₂, y₁₂, 0)^T. Curves plotted in Fig. 6 correspond to r₁₂ = 7.20 () and to r₁₂ = 8.00 (), 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.¹¹ and theoretical observations made by Tasinkevych et al.⁸ 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 cos² θ+ 105 cos⁴ θ (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 r₁₂ = 8.00. This is fully in line with trajectories of colloids measured by Smalyukh et al.¹¹ 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 setup¹¹ 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 ϕ = N_collπr₀²/s_xs_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 n̂₀ [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.¹¹

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 (s_x = s_y = 50 and n̂₀·ê_x = 1). The direction of the far-field director n̂₀ 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̂₀. 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̂(r). However, the Frank free energy is computed only outside the defect core where n̂(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.²⁶ 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.⁴⁰ 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₁₂ and the far-field nematic director n̂₀ if the colloids are sufficiently close to each other.

(2) the position of the minimum shifts monotonically to larger θ if r₁₂ increases.

(3) depending on θ and r₁₂ 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¹¹ 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.

Acknowledgements

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

References

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

---