Colloidal nematostatics

,

Isotropic liquids doped with colloidal particles or just colloids are the classical many-particle systems which, over many decades , have been studied using the methods of molecular and statistical physics, electrolyte theory, physical chemistry, and so on [1,2]. Colloidal particles can be of different physical and chemical origin, but in all cases their interaction via isotropic liquid is of a short range [3], even if the colloids are charged since their electric field is rapidly screened by a numerous couterions. The nematic colloids, that have become widely known just over the last 15 years, are fundamentally different: interaction of colloids mediated by distortions of the nematic director is of a long range. That is why nematic colloids are usually compared not with standard isotropic colloids, but with a system of electric charges. The similarity with the electrostatics has always been a powerful, if not dominating, factor of the development of the field of nematic colloids. But what is the origin of this similarity?
The interaction of particles in a nematic liquid crystal is determined by the elastic free energy (FE) of the director distortions induced by the surfaces of colloids. The director field in the proximity of colloids is governed by the nonlinear Euler-Lagrange equations which are impossible to solve analytically even for a single particle. However, outside these nonlinear near zones, the distortions of the uniform director become weak and take the form of a two-component vector which satisfies the linear Laplace equation. As the scalar electrostatic potential is also described by the Laplace equation, there arises an analogy between the electrostatics and colloidal interaction via the director field. This analogy has been playing a very important heuristic role in the develop-

Introduction
Particles immersed in a nematic liquid crystal interact via the director field n which mediates the distortions induced thereby. The fact that this interaction is of a long range and resembles the electrostatic interaction was first brought to the broad attention by Brochard and de Gennes [8] almost 40 years ago and later by Lopatnikov and Namiot [9]. Over the last 15 years the study of this director-mediated interaction has developed into a new rapidly growing branch of the physics of liquid crystals, i. e., the physics of nematic colloidal systems or nematic emulsions [10]. Depending on the size and nature of colloids, these systems can be very different. The recent development of the field has been inspired mostly by a great diversity of a few-and many-body ordering phenomena in microemulsions with micrometer and submicrometer size colloids [4][5][6][7][11][12][13][14][15][16][17][18]. However, there is also an important class of nematic nanoemulsions with molecular size colloids (dopants) and supramolecular size colloids. One example is a nematic liquid crystal doped with chiral molecules [19] that induce the well-known macroscopic cholesteric ordering [19]. Moreover, nowadays an intensive investigation of thermodynamics [20][21][22][23] and physical chemistry [24] of a nematic liquid crystal doped with a great variety of different solute molecules and micelles [25] is under way. A very special example comprises dye-doped nematic liquid crystals: under the action Mathematically, the similarity between the director-mediated and electrostatic interaction lies in the massless nature of the both theories and the Coulomb-like behavior of the Green functions which derives from it. A massive term in the fundamental functional (e. g., Hamiltonian, FE) has the form of a square of the order parameter describing this system. The theory of electrostatic field ϕ is massless as the term mϕ 2 /2 is absent in the energy functional. It is this property that gives rise to the Coulomb 1/r potential: a finite massive term would have made it short range, i. e., exp(−mr)/r. In turn, the elastic theory of nematic liquid crystals describes the field of the director n which is a unit vector. The massless nature of the director field follows from its definition: the energy term quadratic in n is trivial as n 2 = 1 [19]. As a result, the elastic theory allows the director components that far from the distortion source behave as 1/r and its higher powers which is formally similar to the potential of an electric charge, electric dipole, and so on.
The two systems, however, also have fundamental differences. Electrostatic potential is a scalar described by the linear Laplace (or Poisson) equation. It is the linearity that underlies the definition of the electric charge and its density as the source of electric field. At the same time, n is a vector field described by a linear equation (in the one constant approximation) only in 2 dimensions. Because of this linearity, the deformation source can be straightforwardly established: core of a point defect plays the role of a charge in 2 dimensions [19,[52][53][54]. But in 3 dimensions the field n is described by highly nonlinear equations [19] and point defects cannot be linearly connected with the distortions of n they induce [10,42]. In principle, solution to these nonlinear equations with appropriate boundary conditions determines the director field produced by any source. However, in most cases an analytical solution to this problem is not known and can be found only numerically. Some physically interesting sources, such as a sphere imposing normal director orientation on its surface or defects, are indeed very strong: they produce large magnitude highly nonlinear distortions in their close vicinity that are not tractable in terms of electrostatic analogy [10,42]. If, however, the total topological charge of the source is zero, the deformations are decreasing functions of the distance r from the source so that far enough they become sufficiently weak to allow for a linear description and thus for electrostatic analogy. Then, the theory of a large distance colloidal interaction via the director field, which allows for certain analogy to electrostatic multipole interaction and which we call colloidal nematostatics, consists of two quite different problems. The first problem is to find the form of interaction between different types of colloids that are expectedly elastic analogues of the electric multipoles. The second problem is to express such elastic multipoles in terms of the colloid's shape and anchoring, which implies making a connection between the asymptotic electrostatic-like far director field and its nonlinear source. The fundamental difficulty in solving these problems is that the colloidal nematostatics has no a 33602-3 priori analogue of the Coulomb law. Owing to the Coulomb law, the electrostatic interaction can always be written down in the form of interacting pairs of point-like charges. By contrast, nematic liquid crystal is a field system which has an infinite number of degrees of freedom and, generally speaking, cannot be reduced to a system of interacting pairs. Moreover, in the limit of zero size, a particle is known to produce no director deformations [55]. Hence, nematostatics should consider particles of a finite size.
Since, however, the mathematical treatment of point sources is much easier and more familiar from electrostatics, the theory of interaction in nematic emulsions has been first developed phenomenologically for point-like colloids [10,40,[42][43][44]. The phenomenological energy density, describing a director-mediated 1/r 5 interaction of point-like sources and thus interpreted as the elastic quadrupole-quadrupole interaction, was proposed by Ramaswamy et al. [40]. Later Lubensky et al. [42] further developed this approach to incorporate both the interaction of point-like elastic quadrupoles and the 1/r 3 interaction, interpreted as the interaction of point-like elastic dipoles. The connection of these point-like dipoles and quadrupoles with particular finite-size colloids was established by semi-numerical finding the equivalent multipoles that produce the same far field effect as the real colloids. Namely, the authors of [10,42] showed how to find the equivalent elastic dipole and quadrupole of a sphere with the strong normal boundary condition. The director in the proximity of the sphere is governed by strongly nonlinear equations and is indeed very complicated: the near field contains an accompanying hyperbolic hedgehog or disclination ring [36]. This near field was qualitatively described by an ansatz field based on an expected (visual) similarity with the real director distribution [10,39,42]; from this ansatz, the authors determined the dipolar and quadrupolar asymptotics as functions of the sphere radius and elastic constants.
However, these elastic dipoles and quadrupoles were of a very particular type and did not reflect the general symmetry of nematostatics. Later, the tensorial structure of colloidal interaction was demonstrated by Lev and Tomchuk [43] and Lev et al. [44] directly from the general nematic FE. In this theory, a point-like source was characterized by symmetry of the deformation coat in its close vicinity. It was assumed that the symmetry of the director field in this bulk coating can be represented by distribution of normals to some enclosing surface separating the nonlinear near zone from the linear far zone. The obtained pair interaction potential contained all powers of 1/r : the 1/r term was interpreted as the Coulomb interaction, 1/r 3 as the dipole-dipole interaction, and so on. These potentials had the form of a contraction of certain tensors which, in principle, could be associated with the tensorial characteristics of elastic multipoles.
In particular, the authors of [44] introduced a scalar elastic charge and discussed a Coulomb-like attraction. Its origin, however, remained undetermined. On the one hand, based on the well-known result [19] that the source of the 1/r director behavior is a mechanical torque exerted on the colloid, Lev et al. emphasized that the Coulomb-like attraction can be induced by the vector of external torque Γ [44]. On the other hand, however, it was suggested that, in general, the Coulomb term appears when the elastic coat has neither horizontal nor vertical mirror symmetry. But, if the particle is asymmetric, this is the case without any external torque.
A vector of elastic charge and second rank tensor of elastic dipole were introduced by Fournier as far back as 1993 [56]. For some reason, probably because the multipole interaction was not considered there, this approach remained overlooked 1 and did not receive a further development.
We see that the analogy between the electrostatics and nematostatics has already been greatly developed, but it is not complete when it comes to the 1/r interaction. Indeed, the ultimate source of electrostatic 1/r potential and multipole moments is the electric charge density whose analog has not been established in the nematostatics. In a series of papers [57][58][59][60] we developed the analogy between electrostatics and nematostatics to the level of charge density. In the derived elastic charge density representation, the director components on a sphere, separating the near and far zones, play the role of surface elastic charge density which determines the far field and the elastic multipole moments, section 3. Inasmuch as the director has two independent components, in the elastic density representation there are two densities and two multipole tensors (dyad) in any order. The two-component elastic charge is expressed via the vector of external mechanical torque applied on the particle, and the elastic Coulomb-like coupling between two "charged" particles is found to be proportional to the scalar product of the two external torques. Consequently, the elastic Coulomb interaction can be both attractive and repulsive depending on the relative orientation of the torques: "parallel torques" attract each other whereas "antiparallel torques" repel each other, section 3.6. The developed real space Green function method allowed us to obtain a regular FE expansion in a power series of a parameter (colloid size)/(distance between colloids) and consider the general status of the pairwise approximation to the nematic emulsions, section 3.5. The pairwise approach has its limits: the pairwise interaction potential in a nematic emulsion can be uniquely determined and justified only in the leading order in the small parameter (colloid size)/(distance between colloids) since the next order contains irreducible three-body terms.
Based on the tensorial structure of an elastic dipole and quadrupole, we list all possible types of elastic dipolar particles and important examples of quadrupolar particles, section 4. A connection is discussed between different types of particles and the symmetry of the director field they induce. The proposed classification is illustrated by sketches of some dipolar and quadrupolar colloids. Dipolar colloids are characterized by their isotropic strength, anisotropy, chirality, and some twocomponent vector along the undistorted director (longitudinal component). In particular, we show that a general asymmetric colloid is a chiral elastic dipole. The chirality is eliminated by any symmetry plane, so that quadrupoles cannot be chiral. A particle with two vertical and a horizontal mirror symmetry planes is quadrupolar as all of the four dipolar components vanish. A quadropole can be uniaxial and biaxial. An azimuthal anchoring on the colloid surface can make an otherwise quadrupolar particle a dipole if its easy direction field on the colloid surface has a helicoidal component. In section 5, the theory is applied to the colloid-surface interaction by developing the mirror image method of nematostatics. The method is applied to the charge-wall and dipole-wall interaction.

Linearized director distortions and their source
Consider a 3-d director field n(r), uniform and parallel to the z-axis at infinity, n ∞ = (0, 0, 1), and distorted in a finite number N of particle-like (compact) domains. We refer to such a deformation domain as particle though the distortion therein can be induced by surface of a real particle, by topological defects with zero total topological charge [10,42], or by an external field dying out outside the domain area. Generally, the deformations at the domain center are strong and satisfy nonlinear equations, but sufficiently far they become weak, figure 1. Here the perturbation n ⊥ to n ∞ is small, |n ⊥ | 1, pure transverse, n ⊥ = (n x , n y , 0), and satisfies linear equations. Enclose

33602-5
the i-th domain by a spherical surface S i of radius a i such that outside it |n ⊥ | 1. At the same time, the sphere can be small as compared to the interparticle distances, and the radius-to-distance ratio is a small parameter. The spheres divide the space into the area V in inside and the area V out outside them. In V in the theory is nonlinear, intractable in standard terms (maybe experimentally inaccessible), and the state of i-th particle will be represented by the director distribution n ⊥ (s) on the enclosing sphere S i , s ∈ S i ; a priori this n ⊥ (s) on each sphere depends on the positions and orientations of all particles in the system. The state of the total system will be represented by the director distribution on the total surface B N = N i=1 S i of all N spheres (N -sphere boundary). This is well justified both physically and mathematically. As discussed above, the physical point of view is that the information about the particle that can be considered known is the director distribution close to its surface. The mathematical justification, given below, is that the above boundary condition determines a unique solution to the Laplace equation.
A sphere has the following unique advantages: firstly, it does not introduce any symmetry element so that the symmetry associated with the i-th source is fully determined by the distribution of n ⊥ on the i-th sphere; secondly, it enables us to calculate Green functions which seems to be a task impossible for other surfaces. The final results do not depend on the choice of the surface. In particular, the formulas of this section are valid if i-th particle is enclosed by an arbitrary smooth surface S i .
It is implied in what follows that index t takes values x or y, Greek indices run over x, y and z, index i stands for the particle number, summation over the repeating indices is performed, and vectors s and s refer to a point on the surface S i . For simplicity, we omit the divergence FE terms and work in the one constant approximation [19]. Then the distortion FE functional for the area V out outside all the enclosing spheres is given by The transverse director component n t satisfies the Euler-Lagrange equations which for functional (1) coincides with the Laplace equation n t = 0. Here we assume that the perturbation n t is vanishing at infinity (a sample surface at finite distance will be considered later on). Integrating (1) by parts in the context of this assumption and the equation n t = 0, the FE in V out reduces to the integral over the total surface B N of all spheres: where ν = ν(s) is the unit outer normal to the surface B N at s, figure 1. It is known [65] that solution n t (r) to the Laplace equation in the region V out with the inner boundary B N , that vanishes at infinity, is uniquely determined by the director distribution n t on B N . The formula which expresses the director in V out via the director on the surface reads [65] n t (r) = − where G(r, r ) is the Green function of the Laplace equation [65] which satisfies the following boundary condition: The Green function is the sum G(r, r ) = (4π) −1 |r − r | −1 + I(r, r ) where I(r, r ) is a nonsingular function found from the Laplace equation with the boundary condition (4) 2 . We see that, in the linearized theory, the surface director component n t (s) plays the role of a source of the director component n t (r) in the outer area. While the field equation (

33602-6
solution I(r, r ) to this boundary condition does depend on S i and greatly simplifies for spheres S i . By substituting equation (3) in equation (2) the FE in V out acquires the following important form: where These two formulas form the starting point of our approach. They describe an interaction of the surface elements n t (s)d 2 s and n t (s )d 2 s with the potential U N (s, s ). This potential, however, depends on coordinates of all particles by virtue of the boundary condition on G N which vanishes on all the spheres, not just on the two spheres represented by s and s . Therefore, generally speaking, U N describes a many-particle interaction in accordance with the fact that a field system does not reduce to a finite number of particle-like sources with a pairwise interaction. In what follows we obtain the function I(r, r ) for different numbers of particles, show that n x (s) and n y (s) determine surface densities of what can be considered as a two-component elastic charge, and establish the conditions when formula (6) indeed describes a pairwise interaction of particles imbedded in a nematic medium.

Incorporation of real surfaces and finite anchoring
Here, for a moment, we make a diversion from our main line to demonstrate how the anchoring effect on real surfaces can be incorporated in our approach. Consideration of finite anchoring effects in colloidal nematostatics requires incorporation of a finite-size source, and theories, based on the approximation of point-like sources, encounter difficulties. For instance, the treatment of finite anchoring in [61,62], based on the approach [43,44,63], required an introduction of an artificial director field inside the colloids' material. In our approach, the anchoring effects can be easily included.
Above, for simplicity, the sample surface at finite separations from the colloids, as well as a finite anchoring on real surfaces were disregarded. The above formulas (2)-(5), however, provide a natural basis for straightforward incorporation of these surface effects. Indeed, these formulas reduce the equilibrium elastic FE to the surface integral. But so is the anchoring energy, which thus should be just added to the elastic FE.
If the real surface of the i-th particle is S i and the sample surface is S 0 , the total anchoring energy F a has the form where f i,a and f 0,a are the correspondent anchoring energy densities. Assuming that the director field everywhere (even just outside S i ) is described by the Laplace equation, the sum of the bulk FE (2) and the anchoring energy (7) takes the form of the surface integral which is a functional of the director field on the surfaces S i and S 0 , F N = F N {director field on all of the real surfaces}. This functional is minimized by the standard boundary conditions to the bulk Laplace equation: Substituting here the solution (3) of the Laplace equation and using the definition (6) of the potential U N , one finally obtains following equations for the director on the i-th colloid surface S i 33602-7 and on the sample surface S 0 : The Green function which determines the potential U N via formula (6) now satisfies the zero boundary condition on all of the real surfaces, i. e., Equations (9) are self-consistent integral equations showing that the perturbation of the field at each point of any surface is superposition of the perturbations from all points of all the surfaces. These equations are exact provided the bulk director can be determined from the Laplace equation. Note that system (9) considerably simplifies if the director is fixed on some S k of S i . In this case the equations for the surfaces S k are trivial and disappear from the system (9), while the contributions from S k to other surfaces reduce to the fields of correspondent (known) elastic multipoles.
Thus, the problem of colloidal interaction in all cases, be it with fixed or soft boundary conditions, is reduced to finding the appropriate Green function of the Laplace equation. We are again to stress that the Green function method is the most direct and mathematically most justified and developed way to address the problems of colloidal nematostatics. And, actually, it is this way that makes the essence of the electrostatic analogy (in particular, Jackson's "Classical Electrodynamics" [64] can be as useful as it has been in solving electrostatic problems).
Actual experiments with colloids in a nematic liquid crystal are performed in confined geometries where the bounding surface gives rise to important or even dominating effects such as exponential screening of the 1/r 5 interaction at distances larger than the cell thickness [18]. As in electrostatics, these effects can be addressed by introducing image-multipoles behind the confining surface [59][60][61][62][63]. The image method for finite-size colloids [59,60] is presented in section 3.4. Till then we consider colloids far from confining surfaces.

Elastic charge density and multipoles of a single particle
Consider first a single particle enclosed by a sphere S i of radius a i and set the coordinate origin at its center; all other particles are assumed to be very far. In this case the director field is the unperturbed equilibrium field n (i) t (r) of a single i-th particle which is indicated by the subscript i. The one particle Green function G (i) 1 subject to the boundary condition G (i) 1 (r, s) = 0 for s ∈ S i is known exactly [65,66]: where i(r) = (a i /r) 2 r is location of the image of the point r in the spherical surface which has to be taken with the coefficient a i /r, figure 1. The normal surface derivative of the Green function is obtained in the form [57,65] (ν(s)·∇ s )G Substituting (12) in (3) gives the outer Poisson integral which is exact formula for the director induced outside the single sphere by the director components on its surface, i. e.,

33602-8
Let us introduce the surface density of the transverse director field Then expansion of (13) in inverse powers of the distance r from the origin yields where q t,αβ r α r β ). The form of (14)-(18) prompts the following interpretation. q x and σ (i) y are two separate sources, the dyad, they determine not only the x and y director components outside the single i-th particle, (13), but also two independent tensors, the tensor dyad, for each multipole moment, i. e., q x and Q (i) y , and so on. These multipoles are uniquely determined quantities that characterize the director in the far zone of a single i-th particle unperturbed by other particles. Along with the director n (i) t (r) and its large r asymptotics, given by the r.h.s. of (15), the multipole moments do not depend on the choice of the enclosing surface. Our choice of a sphere surface is just one possibility motivated by its simplicity: as indicated above, it enables one to obtain the formulas for multipole moments in the explicit form given in this section. In particular, below it enables us to express the elastic charge and Coulomb-like interaction solely in terms of the mechanical torques exerted on the colloids.

Torque balance, Gauss' theorem, and elastic charge in 3 dimensions
The fundamental physical quantity of electric charge is purely phenomenological and should be postulated in the theory of elementary particles. By contrast, nematostatics of the director field n allows for introduction of two different charges. Owing to the linearity of the one constant approximation in 2d, the deformation source can be straightforwardly established: core of a point defect plays the role of a charge in 2d [19,[52][53][54]. The integral, expressing the topological invariant, is independent of the integration contour. This property plays the role analogous to Gauss' theorem in electrostatics while the invariant itself plays the role of a conserved charge. As a result, the 2d nematostatics is similar to the 2d electrostatics with its logarithmic potential: disclinations with topological charges of the same sign repel one another and those with topological charges of opposite sign attract one another.
In 3d, however, the analogy between topological defects and charge is completely lost. In 3d, the field n is described by highly nonlinear equations [19] so that point defects, though remain topological invariants, cannot be linearly connected with the distortions of n they induce [10,42].
Here the deformation source is the director distribution in a domain of the size ∼ a [57,58], figure 1. We seek the elastic analog of charge, following de Gennes' idea outlined in [19].
A colloidal system can be subject to external electromagnetic fields. If colloids have a permanent dipole (electric or magnetic) or finite polarizability, these fields exert a mechanical torque on  colloids. Due to its elasticity, a nematic liquid crystal transfers mechanical torque. A torque, exerted on a particle, induces director deformations in the ambient medium. Nonzero deformations result in a nontrivial torque balance in the director field. It is this torque balance that is described by the standard Euler-Lagrange equations minimizing the FE functional. The torque balance implies two torques at any spatial point. Namely, the balance at a given point shows that an elastic torque is applied at both sides of any virtual surface passing through this point, and that these two torques have equal magnitudes and opposite signs [19]. Thus, a torque on a particle, located within a closed surface S, is transferred from inside S to outside S, the total torque being conserved. It means that certain torque density, integrated over an arbitrary closed surface, is equal to the total torque applied inside the surface irrespective of its form. In the static equilibrium, only a transverse external torque Γ ⊥ = (Γ x , Γ y , 0) can be exerted on a particle; a longitudinal torque Γ z cannot be balanced as a rotation about the z-axis does not change the elastic energy, Γ z = 0. In the one-constant approximation, the integral, which expresses the conservation and transfer of a torque Γ ⊥ , is of the form [19] where K is the elastic constant, ε αtρ is the absolute antisymmetric tensor, all indices but t run over x, y, z, and summation over the repeated indices is implied. The integral in the r.h.s. of (19) does not depend on the choice of enclosing surface S, and the equality (19) reminds us Gauss' theorem with Γ t in place of the electric charge. To further justify this connection one notices that integral (19) over a remote surface S vanishes for any term in the expansion (15) except for the first one. Substituting n t = q t /r in (19) and integrating over a sphere gives Γ y = −4πKq x , Γ x = 4πKq y , or The above consideration clearly shows that the Coulomb director field dying out as 1/r can be induced by a colloid if and only if there is an external mechanical torque upon this colloid. In particular, a colloid that has no symmetry planes does not induce the 1/r field if there is no external torque (such a colloid is a chiral elastic dipole whose field is dying out as 1/r 2 , see section 4.2). Thus, the tentative conclusion is that, in 3d, the role of Gauss' theorem and a charge is played, respectively, by the conservation of an elastic torque transferred via the director field, and by two transverse components of an external torque exerted on a particle [59,60]. This is fully justified by calculating the Coulomb-like interaction in the elastic charge density representation [57,58].

The two particle Green function, point-to-point potential and its expansion in a power series of 1/r
Now consider two particles. In contrast to the single sphere case, the Green function G 2 (r, r ) for two spheres cannot be found in the form of a sum of the Coulomb source and a finite number of its point-like images. In [57] we obtained the lowest order terms of expansion of the Green function in a power series of the small radius-to-distance ratio by considering a few successive images of the Coulomb source in two spheres. Here we briefly outline the method and give the main results. Figure 2 shows the two spheres S 1 and S 2 of radii a 1 and a 2 , centered at points o 1 and o 2 and representing two particles. The radius vectors of the points o 1 , o 2 , r, s 1 . . . are denoted by the correspondent bold symbols. The Green function G 2 (r, r ) to the Laplace equation, subject to the boundary condition G 2 (r, r ) = 0 for any r ∈ B 2 = S 1 S 2 , is of the form where I(r, r ) is a field of unknown "image" source that should compensate the field of the Coulomb source 1/|r − r | on the surface B 2 [65]. The vectors i 1 and i 2 represent the first-order images of r, respectively, in S 1 and S 2 ; the vectors i 12 and i 21 are the second-order images of the point r: i 12 is the image of the point i 1 in sphere S 2 , and i 21 is the image of the point i 2 in sphere S 1 . The function I is the sum of all of the image sources. The reflection process can be continued, but we do not need higher images as they result in terms of higher order in the small parameter a/R than the leading order. For r close to the surface S 2 , up to terms∼ O[(a/R) 2 ], the Green function can be written as As seen from general formula (5), the pair interaction potential is the sum U 2 (s 1 , s 2 )+U 2 (s 2 , s 1 ) where U 2 is determined by (6) for N = 2. Finally, the pairwise interaction potential is calculated from the Green function as where, for later convenience, we introduced the factor a 2 1 a 2 2 . Calculation of (23) with the Green function (22) gives the potential in the following form [57,58]: where R = o 1 − o 2 . Now we introduce the unit vector u = R/R in the separation direction, the vectors (directed inward the spheres) of unit outer normals ν 1 (s 1 ) and ν 2 (s 2 ), respectively, to the surface S 1 at point s 1 and to S 2 at point s 2 , and expand (24) in a power series of the inverse separation 1/R. The result can be converted to the following important formula for the pairwise point-to-point potential [57,58]: The four terms shown in this expansion give rise, respectively, to the Coulomb, charge-dipole, dipole-dipole, and quadrupole-quadrupole colloid-colloid interaction potentials which will be derived below.
The arguments s 1 and s 2 of the potential (25) belong to different surfaces. For completeness, we give here the point-to-point potential for two points s 1 and s 2 lying on the same sphere S with the radius a. This potential enters equations (9) and the colloid self-energy is introduced below. This potential is equal to −a 4 lim r−→s (ν(r) · ∇ r )(ν(s) · ∇ s )G 1 (r, s) where G 1 is given in equation (11). The result is 33602-11

Interaction of two particles at large separations, the pairwise approximation and its limits
Now consider the energy of two particles of the previous section, figure 2. The director distribution on S i , i = 1, 2, is now perturbed by another particle. In terms of the elastic charge density (14) this means that the unperturbed single particle charge density σ t (s). Then the total elastic FE F 1&2 of the two particles, obtained from the general many-particle FE (5) for N = 2, is the sum of the following terms: The first three terms are independent of the perturbation δσ. The first term is the self-energy of the unperturbed particle 1 which does not depend on the state of particle 2. Similarly, the second term does not depend on particle 1, so that these two terms do not contribute to interparticle interaction. The third term in (26) is the pairwise interaction FE given by the formula where U 2 is defined by (25). The last term in the FE (26) is the sum of the following corrections. ∆F 12 is the correction to the self-energy F 1 because of the presence of particle 2, and ∆F 21 is the correction to F 2 because of the presence of particle 1. Though the sum ∆F 12 + ∆F 21 appears as correction to the self-FE, it represents interaction as it depends on the locations of both particles. The term ∆ 12 {δσ t } is the total FE contribution due to the final perturbation δσ (i) of the single-particle fields of both particles. We also need the correction due to nonlinearity of the director field. This term, which has the form [42] is not present in the FE (5) of linearized field and should be added individually. The above F 1&2 (26) totally disregards nonlinearity of the exact equations of nematostatics neglecting terms ∼ n 2 ⊥ = n 2 x + n 2 y . In this concern, the natural question arises as to how small should n t be on the integration spheres S i in the formulas (14)-(18) for the multipoles that F an be indeed sufficiently small.
We have shown in [57] that in the system of similar multipoles (only charges, only dipoles, or only quadrupoles), the term F 12 (27) is the principal contribution while all the correction terms in the sum (28) are higher by an order in the small parameter a/R, i. e., The anharmonic correction F an can be neglected as compared to the leading term F 12 if F an (a/R)F 12 , which gives n ⊥ a/R, where n ⊥ is the average magnitude of the transverse director on the integration sphere: n ⊥ ∼ n t n t . This is a rather weak restriction, and n t on the integration sphere in equations (14)-(16) can be quite appreciable. The inequalities (30) and (31) show that the leading order pair interaction potential is F 12 given in equation (27). Moreover, it was shown in [57] that the interaction between colloids via the nematic director field is pairwise only in the leading order because already in the next order there appear irreducible three-particle interaction terms. Now the following important conclusions can be drawn as to the restrictions of the pairwise approach in the case of particles bearing similar elastic multipoles. First, if the particle size and interparticle distance are of the same scale, the omitted terms make the pairwise approximation inapplicable, and one is left with the original field description alone. Therefore, the ratio (a/R) should be small. Second, the pairwise approach is valid with the accuracy up to terms ∼ (a/R)F 12 . Third, the perturbations of the individual multipole moments, induced by other particles and described by the functional ∆{δσ} ∼ (a/R)F 12 , can be adopted within the pairwise approach only in a system of two particles, because only in this case the three-body terms do not appear in this order. This means that the pairwise term F 12 describes the interaction between, respectively, the unperturbed elastic charges, dipoles, and quadrupoles, which are fully determined by the unperturbed director distributions of individual particles. Since, as explained above, these multipole moments are uniquely determined characteristics of an individual particle, the pairwise interaction is also uniquely determined by equations (27). In particular, it is independent of our choice of the enclosing surfaces.
It is known that correlation between two particles via the third one statistically occurs in the second virial coefficient of the thermodynamic virial FE expansion. In our case, the correlations are of static origin and are much stronger: they give rise to FE corrections in the order a/R whereas the statistical mechanical correlations appear in the order (a/R) 3 (the small parameter of thermodynamics is the packing fraction ∝ (a/R) 3 ).

Interaction of two similar elastic multipoles: the Coulomb, dipole-dipole, and quadrupole-quadrupole potentials
Here we consider the leading interaction term F 12 which, at the same time, is the only pairwise term, by substituting expansion (25) in (27). Consider contribution of the first term in (25) to F 12 . Using the definition of elastic charge q t (16) we obtain [57,58] which is obviously the Coulomb interaction. We see that, in contrast to the electrostatic Coulomb law and 2d nematostatics, charge components with similar signs attract while those with different signs repel each other. The relation (20) allows us to express the Coulomb term via the transverse components Γ ⊥ = (Γ x , Γ y ) of the torques exerted upon the particles: This formula shows that the Coulomb attraction is induced when (Γ ⊥ ) < 0, the torques induce a Coulomb repulsion: parallel torques attract while antiparallel torques repel each other. An important conclusion from the formula (33) is that the elastic Coulomb-like interaction can be induced only by torques upon the interacting particles and does not directly depend on their form and anchoring [57,58]. In particular, breaking the horizontal and vertical mirror symmetry of the director field, induced by asymmetric particles, is insufficient for the Coulomb interaction. Calculation of an external torque requires much less information about the source than that contained in the distribution of anchoring strength and easy axes over the colloid surface. For instance, if the anchoring is very strong, then the torque that can be exerted on the colloid by a given (electric or magnetic) field, is determined by the vector of permanent dipole (electric or magnetic) or the polarizability tensor of the particle, and no further information about the anchoring is needed. If the anchoring is weak, then it determines the maximum value of the applied torque.
If the external torques are not applied, the interaction F 12 can be expressed solely in terms of multipoles of the particles. Using definition (17) of the elastic dipole, the potential (27) calculated with the third term of the expansion (25) is obtained in the following form: which can naturally be interpreted as the elastic dipole-dipole interaction. Using the definition (18) of the elastic quadrupole, the potential (27) calculated with the fourth term of the expansion (25) gives the elastic quadrupole-quadrupole potential: Thus, we derived the pairwise interaction F 12 for two interacting charges, dipoles, and quadrupoles. The natural relation between the definitions (14)-(18) and the interaction potentials (32)- (35) supports our idea of the elastic charge density.
Below we will consider the geometrical structure of different elastic multipoles, present a mathematical classification of possible dipoles and quadrupoles, illustrate it by sketching some representative colloid shapes, and consider important examples of the multipole interaction.

Dyads of elastic multipoles and their transformation in the intrinsic (a la isotopic) space
The definition (14)- (18) and interaction potentials (32)- (35) suggest that there are two scalars of elastic charge, two vectors of elastic dipole, two tensors of elastic quadrupoles, i. e., all of the elastic multipoles are dyads. A question arises when this interpretation is correct. Imagine that the director distribution at the source, n x and n y , changes but this change has nothing to do with a rotation about the z-axis, e. g., colloids change alignment at their surfaces dynamically. Then, the new components after the change, n x and n y , have no connection with their old values, n x and n y , via a tensorial transform of the local (i. e., individual for each colloid) reference frame in the real x, y-space. Similarly, the new values of d t,α (Q t,αβ ) have no connection with their old values via a local tensorial transform in the real x, y-space. Such a change can be considered as a local transform in some intrinsic space with the coordinates n x and n y or d x and d y (Q x,αβ and Q y,αβ ), which is individual for each colloid and reminiscent of the isotopic space of nuclear physics. Thus, under such an "intrinsic" transformation, n x and n y do not transform as components of a spatial vector. Then q t is not a spatial vector but rather can be viewed as two independent charges, a dyad; similarly, d t,α (Q t,αβ ) is not a spatial tensor but a dyad of two independent vectors d x,α and d y,α (tensors Q x,αβ and Q y,αβ ). Another obstacle in viewing the objects d t,α and Q t,αβ as tensors is that t runs over 1 and 2 whereas α and β run over 1, 2, and 3. Nevertheless we need to rotate an elastic multipole about the z-axis in the real space to a new reference frame. To this end, in the next section we make it a tensor in the dimension 2 + 1 = 3.

Elastic dipoles in the 2+1 and 2 dimensions
Since the rotation about a homogeneous director n ∞ , which is along the z-axis, does not alter the FE of a particle, all the elastic multipoles should be determined up to an arbitrary azimuthal angular variable φ. Let us consider the subscript α = t with the values 1 and 2 and with the value α = 3 separately. The indices taking values 1 and 2 will be denoted by t or t . Then q t and d t,3 are transformed as a 2d vector, and d tt is transformed as a 2nd rank tensor.
The 2d vector (16) of elastic charge can be formally imbedded into the 2 + 1 = 3 dimensional space by defining the 3d elastic charge vector

33602-14
where Γ = (Γ x , Γ y , 0) is the torque vector. Though "vector of charge" does sound unusual, the tensorial structure of a vector = first rank tensor of elastic charge is quite simple and absolutely clear. An elastic dipole is a more complex quantity. Let us imbed d t,β (17) into the 2 + 1 dimensions. To this end, we introduce a second rank 3 × 3 tensor d α,β with the matrix ||d α,β || whose first and second rows coincide, respectively, with the components of d x and d y ,i. e., In turn, the first two rows can be decomposed into the second rank 2 × 2 tensor D and the 2d vector Consider D and d 3 individually and find their general φ-dependent form induced by a rotation of the particle about the z-axis.
A general 2nd rank tensor D is a sum of its symmetric and antisymmetric parts. In the proper reference frame O 0 , D can be reduced to a special form D 0 with the diagonal symmetric part, i. e., The antisymmetric part is invariant with respect to rotations about the z-axis and has the same form and coefficient C in O 0 and in any other reference frame. A symmetric 2nd rank tensor can be decomposed into an isotropic part ∝ δ ij and anisotropic traceless part. In particular, D 0 (41) (i. e., D in the particular reference frame O 0 ) can be written as where d = (D 11 + D 22 )/2, ∆ = (D 11 − D 22 )/2. Now we want to find the form of the tensor D in an arbitrary reference frame O φ that can be obtained from O 0 by rotating it by an angle φ anticlockwise. Such a rotation induces the following transformation of an arbitrary vector r = (r 1 , r 2 , r 3 ): r → r = (r 1 , r 2 , r 3 ), where r α = R αβ r β and the matrix ||R αβ || is of the form The 2nd rank tensor D 0 is transformed to the reference frame O φ as D tt = R ts R t s D 0ss . In the context of (42) and (43) this finally gives One can easily verify that d 3 transforms as a 2d vector. If d 3 is equal to d 3,y ) in the reference frame O 0 , then d 3,t = R tt d 3,x cos φ + d 3,y cos φ − d 3,x sin φ).

33602-15
The angle φ can also be viewed as a rotation angle of the particle itself in the fixed reference frame O 0 where the dipole has the diagonal form D 0 (44). Then φ is the angle of a clockwise particle rotation from its special diagonal state D(φ = 0) = D 0 to an arbitrary state D: after this rotation the form of the dipole tensor D(φ) will be presented by formula (45). Note that formula (17) provides the following estimates for the magnitude of the coefficients d and ∆. The coefficient d describes an isotropic source, which in our case is uniaxial and produce no spiral component (no chirality): d ∼ a 2 |n t | where |n t | is the average of the absolute value of the transverse director component over the enclosing sphere of radius a. The coefficient ∆ describes an anisotropy related to biaxiality of a nonchiral source: ∆ ∼ a 2 ( |n x | − |n y | ). Since d and ∆ contribute to the transverse components of the dyad (45) normal to the unperturbed director n ∞ , we will refer to them as uniaxial and biaxial transverse dipole strength, respectively. The antisymmetric tensor with the coefficient C describes chiral strength of a source particle. Our idea of chiral source is similar to that by de Gennes [19] and Fournier [56], but it is expressed directly in terms of the director components (see section 4.3.1 below). The vector d 3 is a longitudinal dipole since it consists of the components of the dipole dyad (45) along n ∞ . The formulas (38) and (42) provide a mathematical basis for classifying all the possible types of dipolar colloids. Thus, in general, an elastic dipole can be described by the three coefficients d, ∆, C, and the vector d 3 : the isotropic dipole strength, anisotropy, chirality, and the longitudinal component. Now we discuss a connection between d, ∆, C, d 3 and the symmetry of source particle (which throughout this article implies symmetry of the equilibrium single-particle director field in the vicinity of the particle determined by distribution of the polar and azimuthal anchoring over the particle surface).
Asymmetric particle (general chiral dipole). In the case of a general asymmetric particle, all the three coefficients of the tensor D (44) and the vector d 3 are nonzero. In contrast to the assertion of [44] where the elastic charge was associated with the total absence of mirror reflection planes, such a particle is a general chiral elastic dipole rather than an elastic charge. The C ∞v symmetry (uniaxial particle) (see [67] for symmetry notations). An important particular case of a colloid is a particle with the symmetry C ∞v , i. e., infinite order rotational symmetry about the z-axis plus mirror symmetry with respect to any "vertical" plane passing through the z-axis. Such a particle is uniaxial, i. e., isotropic, ∆ = 0, and nonchiral, C = 0; its dipole is transverse, d 3 = 0. The only nonzero coefficient is d. Here we call such a particle uniaxial; examples of uniaxial particles are sketched in figure 3.
The C ∞ symmetry (chiral helicoid). An azimuthally symmetric particle with the symmetry C ∞ without reflection planes (e. g., a helicoid, figure 4 (a), (b), (c)) is a chiral source with nonzero C, but with ∆ = d 3 = 0. From the physical point of view, the mirror reflection planes of a C ∞v uniaxial particle can be eliminated, and, at the same time, the chirality can be induced, by an azimuthal anchoring with the helicoidal alignment of its easy axes, figure 4 (a), (b), (c). In general, C is eliminated by any single symmetry plane (horizontal or vertical); both d and ∆ are eliminated by a horizontal mirror plane; whereas d 3 is eliminated by two vertical mirror symmetry planes, but can be nonzero for one vertical and one horizontal reflection planes, figure 3 (d). For instance, the well-known bent-core molecules with their planes along the director have a nonzero d 3 . The C 1v symmetry (general biaxial particle), C 2v symmetry (biaxial particle). A C 1v particle with a single vertical reflection plane has no chirality, C = 0, but the vector d 3 is not generally zero. Its tensor D is symmetric and has finite anisotropy ∆ = 0. We refer to such particles as general biaxial. The C 2v symmetry with two mutually perpendicular vertical reflection planes passing through the z-axis eliminates the vector d 3 , d 3 = 0, and particles with this symmetry constitute an important particular case of anisotropic colloids which we call (just) biaxial, figure 3 (b). A uniaxial particle is circular in its horizontal cross-sections. A general biaxial particle is biaxial in the horizontal cross-sections, its director distribution in the x, y planes has mirror symmetry with respect to one of the two principal axes, e. g., it is an isosceles triangle. A just biaxial particle is more symmetric: its director distribution in the x, y-cross-sections has mirror symmetry with respect to both principal axes, e. g., the field in x, y-cross-sections is elliptic. Thus, under rotation of a biaxial particle about the z-axis its elastic dipole transforms as a 2nd rank 2 × 2 tensor (44) with the isotropic fraction of magnitude d, anisotropic fraction of magnitude ∆, and C = 0, in which φ is the angle of the anticlockwise rotation of the particle from the diagonal dipolar state. Very recently we have described possible symmetry operations and point groups of elastic dipoles. These unpublished results can be summarized as follows. There are four pure types characterized by a single nonzero coefficient, the four dipolar singlets, and eight mixed types: single dipolar quintet, single quartet, two triplets, and four doublets. The idea of a single electrostatic dipole is thus left far behind: the nematostatic dipoles are much more numerous and versatile.
A dipolar particle cannot have three mutually perpendicular planes of symmetry as these would make it an elastic quadrupole.

Elastic quadrupoles in the 2+1 dimensions
It turns out that if the dipolar dyad of a colloid is zero its elastic quadrupole Q t,αβ (18) is a much less complex quantity than an elastic dipole. The reason is that a quadrupole particle is highly symmetric. For instance, two vertical and a horizontal reflection symmetry planes eliminate all the four dipolar coefficients and make the particle a quadrupole. Consider such a particle. In the proper reference frame O 0 with the x-and y axes lying in the vertical symmetry planes, the 33602-17 dyad of quadrupole tensors has the following form: To imbed this dyad into the 2+1=3 dimensional space, one can take ||Q 3,αβ || = 0 (similarly to the zero third line in the matrix (37)). Each tensor above is symmetric by the definition of Q t,αβ and depends on a single number, Q x or Q y . If the particle is axially symmetric, then Q x = Q y = Q and we call it a uniaxial quadrupole. The well-known examples of such particles are the Saturn ring quadrupole [36,39,42], induced by a sphere with normal boundary condition, and an onionlike bipolar quadrupole, induced by a sphere with tangential surface director field aligned along the meridians of the sphere, figure 4 (d). If Q x = Q y , particle is a biaxial quadrupole. A biaxial quadrupole is anisotropic in the horizontal cross-sections and can be characterized by the isotropic coefficient Q = (Q x + Q y )/2 and anisotropy ∆ = (Q x − Q y )/2. Choosing the x-axis of the proper reference frame O 0 along the longer symmetry axis of the director near field (which is equivalent to |Q x | |Q y |), we can make both coefficients of the same sign, Q∆ > 0. For instance, an anisotropic modification of the Saturn-ring and onion quadrupoles would be obtained by flattening in the ydirection. Then, for anisotropic Saturn ring quadrupole, Q > 0 and ∆ > 0; whereas for anisotropic onion quadrupole, Q < 0 and ∆ < 0. Obviously, ∆ = 0 for a uniaxial quadrupole. Under anticlockwise rotation (43) of the reference frame by the angle φ from O 0 the components (46) transform as the third rank tensors with Q 3,α,β = 0, i. e., Q t,αβ = R tt R αα R ββ Q t,α β . Performing this calculation gives the dyad (46) in the reference frame O φ : For φ = 0, formulas (47) go over into (46). For a uniaxial quadrupole with ∆ = 0, the tensors Q t,α,β are invariant and (47) remain of the form (46) in any reference frame. The azimuthal anchoring plays an exclusive role for elastic quadrupoles. Indeed, the reflection planes can be eliminated by a helicoidal component of the (tangential) easy axes of azimuthal anchoring, which can make the otherwise quadrupolar particle an elastic dipole. For instance, an onion quadrupole is indeed a quadrupole if the easy directions on its surface are along the meridians connecting the upper and lower poles (boojums). If, however, the field of tangential easy axes is twisted and has a helicoidal component, the onion-like particle becomes a chiral dipole, 4c. This dipole has all the coefficients zero but C which illustrates that C is indeed the chirality. Another interesting example is a particle made of two cones of the same dimensions but with opposite sense of chirality, whose bases are attached to each other along the horizontal symmetry plane. Such a particle is not chiral (due to the symmetry plane), but rather an elastic dipole with nonzero longitudinal component d 3 .
To conclude this section, we note that, depending on its anchoring, a particle of the same shape can be a dipole or quadrupole. For instance, an equilateral triangle-shaped platelett is a dipole if its plane is parallel to n ∞ (planar anchoring), but it is a quadrupole if the plane is normal to n ∞ (normal anchoring).

The electrostatic analogy and tensorial multipole structure of [56]
A vector of elastic charge similar to our q t and the 3×3 tensor of dipole moment were introduced by Fournier in [56]. His analysis of the elastic dipole [56] and our analysis presented above share certain common ideas, too. However, the difference is substantial. Using the techniques proposed by de Gennes, Fournier also introduces the chirality, the isotropic and anisotropic terms, but 33602-18 rather for the rotation vector than for the director components themselves. The two approaches also differ in the choice of the analogy between electrostatics and nematostatics. In our approach the transverse director components n x and n y are considered to be an analogue of electrostatic potential ϕ [39,40,43,44,57,58]. The fundamental reason is that n t enters the FE functional the same way ϕ enters the electrostatic energy. Then, similar to ϕ, which is the canonical variable for the electrostatic field theory, the vector (n x , n y ) is the canonical variable for the nematic field theory. Another way, chosen in [56], is to seek analogy between the director and the electric field, i. e., between n t and the gradients of ϕ. The analogy can be established between the director field induced by an elastic dipole and the electric field of an electric charge (both scale as 1/r 2 ). Of course, technically this can be justified, but such an analogy is no longer direct which results in certain mismatch. In particular, the elastic charge has no electrostatic counterpart (electric charge is already connected with the elastic dipole). We believe that the direct analogy [39,40,43,44,57,58] between n t and ϕ is more natural, deep, and far going than the indirect analogy chosen in [56] which becomes obvious when the multipole interaction is considered.

Particular examples of the elastic dipole-dipole and quadrupole-quadrupole interaction
Now we illustrate the general interaction potential by important examples. First we substitute particular dipole dyads (45) to the general dipole-dipole interaction potential (34) to obtain its particular form for some dipolar particles. To this end, we choose a separation vector R to lie in the xz-plane as shown in figure 5, introduce the polar angle θ made by the unit separation vector u and the z-axis, and the azimuthal angle φ i , i = 1, 2, made by the actual x-axis and the x-axis of the proper reference frame O (i) 0 of the i-th particle, figure 5. Then, the separation vector lies in the xz plane, u = (sin θ, 0, cos θ), and the dipole diad of the i-th particle is given by formula (45) in which all the parameters are indicated by index i. The i-th C ∞v uniaxial dipole is characterized by ∆ = 0, d (0, d i , 0). The interaction potential of two uniaxial particles is of the form which, in particular, describes the interaction of topological dipoles [42]. The i-th C ∞ uniaxial chiral dipole is characterized by ∆ = 0, d , d i , 0). The interaction potential of two uniaxial chiral particles is of the form

33602-19
This formula is incapable of demonstrating the full effect of chirality which can be described only if the restriction of uniform unperturbed director is removed. Then the effect of chirality is expected to be a spiral cholesteric-like rotation of the director field. The dipole dyad of an i-th biaxial dipole is obtained from (45) for cos 2φ i , 0). The interaction potential of two biaxial dipoles has the form Now we substitute quadrupole diad (47) to the general quadrupole-quadrupole interaction potential (35) to obtain its particular form for uniaxial and biaxial quadrupolar particles. The i-th quadrupole is determined by coefficients Q i and ∆ i , (47). Two biaxial quadrupoles interact via the potential The potential for two uniaxial quadrupoles is obtained by setting ∆ i = 0, which gives This formula differs from that derived for Saturn ring quadrupoles in [42] only by coefficient. General formulas for the interaction potentials (32)- (35) and their particular forms derived in this Subsection solve the first problem of colloidal nematostatics formulated in the Introduction: they provide the form of interaction between different types of colloids. However, these formulas can be used if the material parameters of the multipoles are known. The problem of finding these parameters in terms of the shape of the colloid and anchoring is entirely different from the first one and can be addressed only numerically. Nevertheless, formulas (17) and (18) greatly simplify the calculation of a dipole moment. For instance, the coefficient d for the topological dipole and quadrupole [42] mentioned in this review more than once, was calculated by using ansatz (29) of [42]. Although this anzatz describes the near field only qualitatively, integrating it over the sphere, enclosing the pair, closely reproduces the result which was obtained in [42] by tailoring the fields in the far and near zones [57,58]. Notice, that equation (17) was used in [16] to estimate the dipole moment of liquid droplets trapped at a nematic-air interface.

The mirror image method of colloidal nematostatics
The above formulas can be used to solve boundary problems similar to those of electrostatics. The simplest boundary problem is the interaction of an elastic multipole with a plane surface bounding the nematic sample and imposing a fixed homogeneous director alignment. Here we consider this problem for the elastic charge and dipole. The field of colloidal emulsions deals with macroscopic distances, and here we are interested in macroscopic colloid-wall separations of an order of micrometers. Until recently, the problem of calculating a mechanical torque on a particle near a nematic liquid crystal surface and its energy has been specified by a particular shape, anchoring, and orientation of the particle, and thus could be addressed only numerically [68,69]. The theory [57,58] developed above considerably simplifies the problem. The specific parameters of a particle now enter the problem via its multipole moment, and the problem reduces to the interaction between a wall and an elastic multipole. Now this problem can be solved analytically in a universal form by means of an image method specific to nematostatics. Using the mirror image method of electrostatics as a guiding idea, we developed a mirror image method in nematostatics for arbitrary  director tilt at the wall [59,60]. A wall is shown to induce a repulsive tilt-dependent 1/R 4 force on an elastic dipole and its reorientation to the minimum energy alignment. The external torque, however, induces the elastic charge in this colloid and triggers switching to the 1/R 2 repulsion. The calculations demonstrate that the dyadic nature of an elastic dipole is essential.
Consider a single particle at distance h from a plane surface of a nematic liquid crystal sample (wall). We assume that the anchoring of the wall is strong, the director alignment in the sample far from the particle is homogeneous and parallel to the z-axis, n ∞ = (0, 0, 1), but the angle θ it makes to the surface normal is arbitrary, figure 6; θ = π/2 and θ = 0 correspond to the planar and homeotropic surface director alignment, respectively. To justify the linearized theory, h is assumed to be large compared to the particle size ∼ a. As the director on the wall is fixed, the boundary condition consists of two equations which have to be satisfied for any point r wall of the wall. This boundary condition (53) can be satisfied by placing an image-particle with the multipole moment of the same order on the other side of the wall. Consider an elastic multipole at a distance h from a wall and its image multipole located at the mirror point behind the wall at the distance R = 2h from the real particle, figure 6. Variables, attributed to the particle and its image, will be indicated by index 1 and index 2, respectively. The director field, which determines the multipoles, is fixed at the spheres S 1 and S 2 of radius a enclosing, respectively, the particle and its image. To the order a/R, the transverse director components n t on the wall are given by the sum [58] n t (r) = n (1) where n (1) t is the far field (15) of the single particle, where o 1 and o 2 are the radius vectors of the two centers. Equation (54) explicitly shows nonadditivity of the fields produced by the two particles: n t = n (1) t on S 1 (r 1 = a) and n t = n (2) t on S 2 (r 2 = a); the additivity n t (r) ∼ n (1) t + n (2) t takes place only when both a/r 1 1 and a/r 2 1, i. e., far from the particles. We set the reference frame with the onset on the wall between the particles so that the xz-plane (with z-axis along n ∞ ) is perpendicular to the wall, and the y-axis is normal to the figure plane, r wall = (−z cos θ, y, z sin θ), and r 1 = (−z cos θ + h sin θ, y, z sin θ + h cos θ), r 2 = (−z cos θ − h sin θ, y, z sin θ − h cos θ).
The particle-image separation vector is R = o 2 − o 1 = Ru, where u = (sin θ, 0, cos θ). In the context of (54) and the equality r 1 = r 2 , the condition n t = 0 on the wall yields These equations determine the components of the image multipole. We will find the image charge and image dipole for arbitrary tilt θ and then, by substituting its components in the interaction energy (32) or (34), calculate the interaction between the colloid and the wall. Below we consider a charged, nonchiral C ∞v uniaxial, biaxial, and general dipolar particle individually.

Elastic charge-wall interaction
The charge q (1) t can be induced by an external field exerting the torque Γ with the components Γ y = 4πKq (1) x and Γ x = −4πKq (1) y ; the image-charge q (2) t is on the other side of the wall at the distance h from it, figure 6. The director field of a single elastic charge is n t , t = x, y. Thus, the particle and its image have opposite charges which corresponds to Γ figure 7. Since two opposite elastic charges repel one another, the elastic charge-wall interaction is repulsive. The repulsion force is obtained from the interaction Figure 7. Elastic charge at a wall with fixed planar (left sketch) and homeotropic (right sketch) director alignment. Elastic charge qt induced by an external torque Γ and its image −qt induced by the image-torque −Γ. The director at the wall remains unperturbed and equal to n∞. energy (32) of the torques Γ ⊥ and −Γ ⊥ by differentiating with respect to R at R = 2h, i. e., The result depends on the direction of n ∞ and thus on the surface tilt only via the trivial relation Γ ⊥ = Γ − (Γ · n ∞ ).

Interaction of a uniaxial C ∞v dipole with a wall
Now consider the interaction between a wall and an elastic dipole dyad (d x , d y ), (45). Consider first the simplest and practically important case of a uniaxial C ∞v particle with the isotropic nonchiral dipole dyad:

33602-22
Substituting (58), (59) and (55) in (60) yields the equations which can be reduced to the form This system is solved by d x = (−d cos 2θ, 0, d sin 2θ), Substituting (62) in (48) with u = (sin θ, 0, cos θ) yields Differentiating (63) with respect to R at R = 2h, we obtain the force with which the wall repels a uniaxial dipole: Now consider two important particular cases of the planar and homeotropic wall. For the planar wall θ = π/2, the image has the form d For the homeotropic wall θ = 0, d x = (−d, 0, 0), d y = (0, −d, 0), figure 9, and the repulsive force has the magnitude The force F 1d−wall,hom is 1.5 times weaker than F 1d−wall,planar . In this concern we would like to emphasize that the results, obtained above by means of nematostatics of [57,58], show that even in the simplest case of a homeotropic and planar wall, the dyadic nature of the elastic dipole is essential. Indeed, a single dipole component is insufficient to satisfy the two boundary conditions (53) on the wall. Therefore, phenomenological theories, such as [42], which deal solely with a uniaxial dipolar particle and describe it by a single coefficient d, cannot be used in the problem of a particle-wall interaction. Nevertheless, let us compare our result with what could be "naively" obtained from the phenomenological approach in the simplest cases. When the poles of a single dipole interchange, the dipole changes its orientation to the inverse one and is described by the coefficient (−d). In this situation, the image dipole can be chosen only in the form of a singlecomponent dipole with the same or inverted sign. In the homeotropic case, the image dipole 33602-23  d (2) = −d allows us to eliminate the particle-induced perturbation n t on the wall and thus to satisfy the boundary condition (56). In the planar case, however, the y-equation cannot be satisfied: the choice d (2) = −d is unacceptable as it would result in attraction which is obviously incorrect, whereas the choice d (2) = d satisfies only the x-equation (56). Nevertheless, this last choice is the only possible one in the sense that at least it results in a repulsion from the wall [51]. Then, using (48) with |d (1) d (2) | = d 2 , one obtains for the planar (θ = π/2) and homeotropic (θ = 0) cases, respectively, the following formulas: Apart from the difference in the coefficients, the "naive" prediction is that the repulsion from a wall with homeotropic anchoring is twice stronger than that from a wall with planar anchoring, F 1d−wall,hom / F 1d−wall,planar = 2. This is in contrast to the above exact result F 1d−wall,hom / F 1d−wall,planar = 2/3.

Elastic charge and simple uniaxial elastic dipole in the gravitational field: a sketch of experiment
The above prediction can be verified in an experiment similar to that recently reported by Pishnyak et al. [51] who studied the interaction between colloids of the topological dipole type [42] and a nematic liquid crystal surface with a planar anchoring. In a thick cell with horizontal surfaces, the distance h from the lower surface is set in balance between the attraction due to the weight of the colloid and the elastic dipole-surface repulsion. The ratio h h /h p of the equilibrium distance h h of a topological dipole from a lower surface with a homeotropic anchoring to the distance h p from that with a planar anchoring can be estimated experimentally and compared both with the ratio h h /h p = 4 2/3 0.9 obtained above, and with h h /h p = 4 √ 2 1. 2, following from the "naive" prediction (67). Note that the value of the coefficient d in our theory, which does not enter the ratio h h /h p , is 1/3 of that calculated numerically in [42] (see [57,58]).
Another interesting effect to observe is "charging" of a colloid, i. e., inducing an elastic charge by exerting external torque [e. g., by applying an electric (magnetic) field if the particle is ferroelectric (ferromagnetic)]. A torque on a colloid creates an elastic charge which is repelled from the wall with a force (57). In the case of a dipole-type colloid, an external torque would switch the h −4 dipolewall repulsion to a h −2 charge-wall repulsion which can manifest itself, in particular, in a sharp increase of the equilibrium distance from the wall. Such an effect was observed by Lapointe et al. in 2004 [70]. A ferromagnetic Ni wire with the planar anchoring along the long-axis in zero field was aligned with the director and rested on the lower substrate of a planar nematic cell. A torque exerted by an external magnetic field turned the wire thereby inducing the elastic distortions of the director field. As a result, the wire repelled from the substrate and levitated over it 3 .

Interaction of a general dipole with a wall
The mirror-image method can be applied to dipolar particles of arbitrary shape. In general case, dipole dyads of a real particle and its image have all the three components, d It is important to remember that the components of the dipole dyad entering this formula depend on an arbitrary angle φ, (45). This fact is trivial in a homogeneous space as the particle energy is degenerate in φ, but the presence of a wall breaks the azimuthal symmetry around the z-axis and removes the degeneracy. This means that the angle φ changes to minimize the interaction energy with the wall. In other words, a particle, approaching the wall from a large distance, will turn to assume some particular orientation with φ = φ m . Below we consider an important particular case described by formula (70) with an explicit φ dependence.

Interaction of a biaxial dipole with a wall
From (45), dipole dyad of a biaxial particle is obtained as d Minimization with respect to φ (the free energy and the force have the same φ-dependence (71)) yields: i) φ is arbitrary for a uniaxial C ∞v particle with ∆ = 0, which is obvious since such a particle is invariant with respect to rotation about the z-axis; ii) φ is arbitrary at the homeotropic wall, θ = 0, which is obvious since such a wall does not break the symmetry about the z-axis, Fig.  5; iii) for 0 < θ π/2, φ = 0 if ∆ < 0 and φ = π/2 if ∆ > 0; this means that the longer axis of the x, y-cross-section of the particle turns parallel to the wall which obviously saves the elastic energy. Particles producing director deformations with two symmetry planes passing through the zaxis, referred to here as just biaxial, can be exemplified by objects of the following geometries: the topological dipole of [42] created by an ellipsoid with the normal surface anchoring rather than by a sphere; an elliptical cone with its long axis along n ∞ and planar anchoring at its lateral surface. In all of these cases the equilibrium orientation of the long axis in the transverse x, y-cross-section is parallel to the wall, and thus the repulsion from a wall is accompanied by a reorientation about the homogeneous director. The only exclusion is the case of homeotropic wall where the only effect is repulsion.

Conclusion
The theory developed above deals with two colloids. However, as the director-mediated interaction is of a long range, the many-particle effects can dominate in the large scale phase behavior of nematic colloids. The cholesteric spiral induced by molecular chiral dipoles is the best known example [19,56]. The theory consistent with the many-particle effects should adopt a slow spatial director variation such as the above cholesteric spiral. Such a theory is our goal for the near future, and the two-particle theory presented above is going to be the ground for building this manyparticle generalization. To conclude this review we would like to bring to the readers' attention three known systems in which many-particle effects are or expectedly should be essential.
The first one is a kind of elastic plasma. Dye molecules aggregate in a nematic solvent and form columnar aggregates [27]. Their long axes make a certain nonzero polar angle to the director whereas azimuthal angles are distributed isotropically due to thermal randomization. When an external electric field is applied, the torques exerted on the colloids are also distributed homogeneously in the plane normal to the director: two aggregates, whose azimuthal angles differ by π, experience opposite torques (e. g., one is turned clockwise, another one anticlockwise). In terms of the elastic charge, this is an elastic plasma: the external torque charge the columnar colloids while the temperature makes the number of charges with opposite signs equal, i. e., the total charge is equal to zero. The elastic charges interact via the Coulomb potential (32) which results in a collective macroscopic director reorientation. This effect in [27] was obtained in a phenomenological way since the potential (32) was at the time unknown. The effect is awaiting for its many-body consideration in terms of elastic plasma.
Another system is a nematic doped with ferromagnetic colloidal particles [29]. Under the action of magnetic field such colloids will become elastic charges, too. We expect that the many-body interaction effects are important to the physics of such ferrocolloids. Interaction effects have not been considered so far.
The third system is formed by colloids trapped at a nematic-air interface [14 -16]. The colloids form hexagonal and quasi-hexagonal lattices with different periods which can coexist [15,16] and be transformed into each other under certain conditions [16]. The lattices are stabilized by the elastic repulsion between the dipolar colloids and the capillary attraction induced by a vertical force on the colloids [16]. This force is definitely due to the elasticity of the nematic liquid since stabilization holds only while the liquid is in its nematic phase and disappears in the isotropic phase. The manybody effect in the capillary attraction is shown to be by 4 − 5 orders of magnitude stronger than the pair interaction [71]. At the same time, the many-body effects of the long-range elastic dipolar repulsion has not been considered.
In our opinion, the next step in the theory of nematic colloids is to incorporate the long-range many-body colloidal interaction for the purpose of describing spatially inhomogeneous phases.