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

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 mesoand 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 301. 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.


Introduction
Liquid crystals are fluids made of molecules that lack spherical symmetry. Instead, their molecules contain elongated, rigid cores that form nematic liquid crystals, or disk-like cores that produce discotic liquid crystals or even more complex shapes. This simple property of anisotropy produces a myriad consequences for the optical, mechanical and thermodynamic properties of liquid crystals. For example, as the temperature is decreased they undergo a series of phase transitions where the symmetry of their state is spontaneously broken. Starting from a high temperature isotropic fluid, where all orientations are equivalent, to a nematic state where orientational order is broken, to a smectic phase, where positional order is broken in one dimension.
The molecular anisotropy produces effects on a macroscopic scale. In the nematic state, all molecules tend to align in the same direction, called the fluid's director. 1 But any localized deviation of molecular orientations from the director will cause a restoring, elastic force. However, a single global orientation cannot be satisfied for all boundary conditions. A colloid immersed in the nematic fluid causes a specific preferential alignment of the liquid crystal molecules on its surface, which is termed anchoring. A spherical colloid has a symmetry incompatible with a global nematic director. Thus, defects in the orientational field will arise. These topological defects are points or lines that minimize the free energy of the liquid crystal subject to complex boundary conditions. Depending on details of the fluid and the anchoring of its molecules at the colloid, the orientational field can be of such dazzling complexity that experts are just beginning to unravel its structural details. 2 If multiple colloids are immersed in a nematic liquid crystal, the distortions of the local director field also give rise to effective interactions between the colloids mediated by the nematic host fluid. 3 These interactions may be used to self-assemble the colloids into supramolecular entities in a controlled (i.e., directed) manner. In this way ordered assemblies of colloids of an enormously rich structure with a great variety of symmetries may be built that would not otherwise exist without the ordered nature of the host phase. 4,5 These self-assembled colloidal structures are also of practical importance, as they can be used to produce photonic crystals. 6,7 The forces guiding colloidal self-assembly result from a complex interplay of molecular scale ordering, mesoscopic elastic interactions, and large scale topological arguments. While the framework of the Landau-de Gennes theory is appropriate for mesoscale effects, the elastic Frank-Oseen free energy is appropriate for long-range interactions. 8 However, there is a gap between the microscopic, molecular information and the coarse-grained approaches used at the meso and macroscopic scale. No consistent method exists to transfer physical information from the molecular scale up to the scale where elastic forces and nonlinearities in the nematic interactions occur. Here, we pursue the goal of bridging this gap. The task is not idle, as there are physical situations that have so far eluded a precise explanation at the molecular level. As a test case to demonstrate power and validity of a novel hybrid approach developed in this work we take on the problem of colloidal self-assembly. Poulin and Weitz 9 found experimentally that in a nematic phase the colloidal center-to-center distance vector r 12 forms a ''magic'' angle of y E 301 with the global director n 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 501 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 y . . . 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 y E 301 for certain r 12 = |r 12 |; this minimum shifts to larger y as r 12 increases. However, to date and to the best of our knowledge no molecularscale 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 y. 11 For example, at y E 301 the colloids attract each other whereas at y = 01 and 901 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.

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 0 . The following three sections are devoted to introducing the various constituents of our model and to specifying their interactions with one another.

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 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 (j iso ) and an anisotropic part (j anis ) where the latter depends on the orientations o i and o j of the mesogenic pair. In fact, o i,j = (y i,j , f i,j ) where y i,j and f 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 e mm is the depth of the attractive well, s is the van der Waals diameter of the isotropic core, and j rep and j att represent repulsive and attractive contributions to j iso , respectively. To derive an expression for the anisotropic contribution in eqn (1) we follow Giura and Schoen. 12 From a lengthy but relatively straightforward derivation these authors show that In eqn (3) the anisotropy function is given by where û(o 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 2 (x) = 1 2 (3x 2 À 1) is the second Legendre polynomial and the (dimensionless) anisotropy parameters 2e 1 = Àe 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.

External field
For a suitable choice of thermodynamic state parameters the host phase introduced in Section 2.1 exhibits isotropic and nematic phases. 12 Unfortunately, when a nematic phase is forming in the bulk it is impossible to determine beforehand the direction of n 0 . In fact, there is an infinite number of easy axes 15 with which n 0 may align in the bulk nematic phase.
To predict and control n 0 it has therefore become customary to place the liquid crystal between solid substrates with specially prepared surfaces that tend to align mesogens in their vicinity in a desired way. 16 Because orientational order in a nematic liquid crystal is a long-range phenomenon, substrate-induced alignment of mesogens allows one to control n 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 where Dz i = z i AE s z /2 is the distance of mesogen i from the lower (+) and upper (À) substrate plane, respectively. Hence, j ms may be viewed as a local, orientation dependent external field, the strength of which is controlled by e ms . Throughout this work we take e ms /e mm = 5.00. The value of e ms is chosen for two reasons. First, it guarantees a sufficiently strong alignment of mesogens with the surface so that fluctuation of n 0 over the course of the simulations are negligible. Second, e 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 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 |û(o i )Áê x | o 1. In other words, g(o i ) may be viewed as a mathematical ''device'' mimicking aligning substrates in experimental setups. On account of its definition in eqn (6), g(o i ) allows one to more or less pin n 0 to the x-axis on average where |n 0 Áê x | C 1 due to thermal fluctuations.

Colloidal particles
In addition, colloidal particles are immersed into the nematic liquid crystal. These colloids are spherical in shape, where r 0 = 3.00s 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 18 we take the mesogen-colloid interaction potential to be given by where 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 0 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 e mc which we maintain constant so that e mc /e mm = 3.50. Again, this value has been selected on the basis of the same rationale given in Section 2.2.
In eqn (7), Z is the inverse Debye screening length and parameters a 1 and a 2 have been introduced to guarantee that the minimum of j mc remains at a distance r 0 + s from the colloid's center of mass and to maintain the depth of the attractive well at e mc irrespective of Z. 18 Throughout this work we fix Zs = 0.50.
Last but not least, we introduce another anchoring function in eqn (7) which we take to be given by 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, 15 locally planar fashion with respect to the colloid's surface.

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.

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 0 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 as where S(r) is the nematic order parameter, n a (r) is the a-component of n(r), and d ab 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 where the Landau-de Gennes contribution is given by using Einstein's summation convention. Eqn (11) is essentially a Taylor expansion of the free-energy density in terms of the order parameter tensor at the isotropic-nematic (IN) phase transition. In eqn (11), A, B, and C are unknown expansion coefficients depending only on density r and temperature T.
Assuming uniaxiality, we employ the identities Q ab Q ba = 3 2 S and Q ab Q bg Q ga = 3 4 S 3 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) is a similar LdG free-energy density of the nematic host phase obtained for the same T and r 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 Df 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][20][21] The latter is defined through the inequality S(r) 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 Df 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 freeenergy 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 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 DF. This is in particular so if the structure of S(r) and n(r) in thermodynamic equilibrium is rather complex. 24 Second, there is no way within the framework of the continuum approach to determine the material constants A, B, C, and K (or L). Here one usually resorts to experimental information which is available for a few liquid crystals. 25 In closing, we emphasize that the expressions given in eqn (12) and (14) correspond to the same ground state. For example, in the absence of any perturbation, that is if S(r) = S b and n(r) = n 0 , D f LdG = f el = 0.

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 DF [after an integration of Df (r) over volume, see also Section 4.3].
Thus, to eventually compute DF we begin by introducing the local alignment tensor at the molecular level where 1 is the unit tensor and h. . .i denotes an ensemble average. 18,26 In eqn (17), r(r) is the local number density. Because Q(r) can be represented by a 3 Â 3 matrix it has three eigenvalues which we obtain numerically using Jacobi's method. 27 We take the largest eigenvalue of Q(r) as the local nematic order parameter S(r) and the associated eigenvector as the nematic director field n(r). However, to compute D 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 1 , K 2 , and K 3 associated with bend, twist, and splay deformations of n(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 C K 2 C 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 C 1.66e mm s À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 ).
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 x = cos y, and y is the azimuthal angle if we take the z-axis of our coordinate system to coincide with n 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 mesogenmesogen interactions to the free energy, respectively. 32 The integrand on the right-hand side of eqn (18) contains the orientation distribution function a(x). It depends only on y due to the uniaxial symmetry of the nematic phase and is normalized according to which implies that in the isotropic phase a(x) = 1 2 . Again, because of the uniaxiality of the nematic phase we expand a(x) in terms of Legendre polynomials {P l } via Inserting this expression into eqn (18), expanding the integrand in a Taylor series around x = 0 (i.e., at the IN phase transition), and retaining in this expansion terms up to O (x 4 ) allows us to recast eqn (18) as if we limit ourselves to the leading term l = 2 in the expression for x and use S 2 = S b . In eqn (21), a = 2rk B /5 and T* = À5ru 2 /2k B where the latter is the temperature at which the isotropic phase becomes thermodynamically unstable. Assuming Df n = Df 0 one easily obtains by comparing terms of equal power in S b in eqn (21) with corresponding ones in eqn (13). Hence, B o 0, C 4 0, and A changes sign at T = T* as it is to be expected at the IN phase transition. 1 3.2.3 Elements of finite-size scaling theory. Unfortunately, the expression for T* given in the preceding section depends on 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 = À 32pe 1 e 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.
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 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 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 L Àd where L is the linear extent of a system of dimension d. 34 Moreover, it has been demonstrated by Vollmayer et al. 35 that the ''distance'' of a unique cumulant intersection from the point at which the phase transition would occur in the thermodynamic limit scales as 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.

Numerical details
Our results are based upon MC simulations in a specialized isothermal-isobaric ensemble discussed in detail elsewhere. 18 In this ensemble a thermodynamic state is specified by N, T, the ratio of side lengths of the simulation cell in the x-and y-directions s x /s y , the distance s z between the solid substrates, and the transverse component P J = 1 2 (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. 36 We generate a Markov chain by a conventional Metropolis protocol where it is decided with equal probability whether to displace a mesogen's center of mass by a small amount or to rotate the mesogen around a randomly chosen axis. All mesogens are considered sequentially such that a MC cycle commences once each of the N mesogens has been subjected to either a displacement or rotation attempt.
Our results are typically based upon 1.5 Â 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-ofmass separations r c = 3.00s. 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.50s.
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 e mm , length in units of s, and temperature in units of e 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 r C 0.92 for which the host phase is sufficiently deep in the nematic phase indicated by S b C 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 0 .

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.
Another well-known finite-size effect is visible in the isotropic phase (i.e., for T 4 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 j 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 C 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 r nor T. This is similar to Maier-Saupe theory where a similar universal value is found. 38 However, more recent MMF DFT calculations showed that instead S IN exhibits a temperature dependence such that a limiting value of S IN is approached from above if T IN is sufficiently high. 12 To estimate T IN (and therefore S IN ) we resort to FSS theory and compute g 2 from eqn (25) using S n b ¼ S n b 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 hS b i p N À1/2 and that of hS b 2 i p 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 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 Df 0 from eqn (13) in Fig. 3. As expected, the absolute minimum of Df 0 corresponds to the isotropic phase (S b = 0) for T 4 T IN . Exactly at T = T IN , Df LdG exhibits a double minimum, one at S b = 0, the other one at S b 4 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 4 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 Df 0 exhibits a saddle point at S b = 0 indicating that the isotropic phase is unstable for all T r T*.

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 welldefined boundary. As already mentioned in Section 3.1, special precaution has to be taken to treat the contribution of the defect core to Df (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 E 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. 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 4 3 pR 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 C 1.80 [see Fig. 4(c)] we obtain f core E 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) where V core = {r|S(r) 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 ds = 0.2 and count the number of cubes in which S(r) r S IN . The total volume of all these small cubes is then equal to V core . In a similar fashion we compute D f (r) for small equally sized volume elements and obtain DF through a three-dimensional integration of D 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. 20 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. 20 who show that even for rather disparate n(r) the free energy of the defect core is more or less the same such that F core p V core to a reasonable degree.

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 y = 01 between the intercolloidal center of mass distance r 12 and the far-field nematic director n 0 a toroidal defect structure exists surrounding the point of contact between the colloids. As y increases this torus is first ''ripped apart'' and eventually a handle-like structure forms at y = 901. 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 0 except in the immediate vicinity of the colloidal pair. As one approaches the colloids, n(r) is deformed with respect to n 0 .
Finally, we plot in parts (g)-(i) of Fig. 5 the local biaxiality order parameter z(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 l À = À[S(r) + z(r)]/2 and l 0 = À[S(r) À z(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, z(r) is largest.
This also offers another route to treating 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 y. Hence, it seems intuitive to introduce an associated change in effective free energy DF eff = F el + F core + DF LdG À 2DF B (28) where DF B is the change in free energy associated with a single Boojum defect. In computing DF 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 DF B is roughly of the same order of magnitude in the range of k B T in agreement with plots presented in Fig. 6 and in Fig. 7. Data for DF eff plotted in Fig. 6 exhibit a number of interesting features. First, one notices that depending on y, the sum F el + F core + DF LdG may be larger or smaller than 2DF B and therefore DF 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 DF eff shifts to larger angles y 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 y 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 p 9 À 90 cos 2 y+ 105 cos 4 y from which y min E 491 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 DF eff is repulsive. In particular, such a repulsive region exists for angles y \ 701 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 DF 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 DF 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 DF 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 DF eff between the four nearest-neighbor grid points in a bilinear fashion to get DF 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(ÀDF eff /DF B )].
For two packing fractions f = N coll pr 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 y E 301 with n 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

Discussion and conclusions
By means of a novel hybrid approach we investigate the selfassembly 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 liquidcrystal 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 0 . They depend on the arrangement of the colloids in space.
We are treating both the local lowering of nematic order and the distortion of the director field at mean-field level assuming that the former can be adequately accounted for by LdG theory and the latter within the standard expression for the Frank elastic free energy. The mean-field character is reflected by the fact that we compute LdG expansion parameters a, B, and C from mean-field DFT and that the elastic contribution to the change in free energy does not account for fluctuations of n(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. 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 secondorder 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 y E 301 between the intercolloidal distance vector r 12 and the far-field nematic director n 0 if the colloids are sufficiently close to each other.
(2) the position of the minimum shifts monotonically to larger y if r 12 increases.
(3) depending on y 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, canonicalensemble 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.