Anisotropic electrostatic screening of charged colloids in nematic solvents

See allHide authors and affiliations

Science Advances  27 Jan 2021:
Vol. 7, no. 5, eabd0662
DOI: 10.1126/sciadv.abd0662


The physical behavior of anisotropic charged colloids is determined by their material dielectric anisotropy, affecting colloidal self-assembly, biological function, and even out-of-equilibrium behavior. However, little is known about anisotropic electrostatic screening, which underlies all electrostatic effective interactions in such soft or biological materials. In this work, we demonstrate anisotropic electrostatic screening for charged colloidal particles in a nematic electrolyte. We show that material anisotropy behaves markedly different from particle anisotropy. The electrostatic potential and pair interactions decay with an anisotropic Debye screening length, contrasting the constant screening length for isotropic electrolytes. Charged dumpling-shaped near-spherical colloidal particles in a nematic medium are used as an experimental model system to explore the effects of anisotropic screening, demonstrating competing anisotropic elastic and electrostatic effective pair interactions for colloidal surface charges tunable from neutral to high, yielding particle-separated metastable states. Generally, our work contributes to the understanding of electrostatic screening in nematic anisotropic media.


Highly charged biopolymers like DNA and filamentous actin are just two of many examples of biological relevance of electrostatic interactions that are screened by counterions under physiological conditions. In soft condensed matter, similar effects allow for exploiting electrostatic interactions between particles in defining colloidal self-organized superstructures that they can form and, even more importantly, enabling the very existence of metastable colloidal systems. These structures range from regular—linear, two-dimensional (2D), and 3D crystalline—structures to amorphous structures and have broad relevance, including in photonics, optics, paint, and food industry. Screened electrostatic interactions have been studied quite extensively in systems where the solvent is isotropic. Findings include colloid stabilization by charge (1, 2), the description of short-range liquid order in scattering experiments (3), measurements and calculations of pair interactions between two charged-screened particles (47), ion transport (8), the influence of external fields (911), and predicting phase behavior, such as demixing (12, 13) and crystallization (14, 15). Furthermore, there are many theoretical studies on many-body effects (1618), charge regulation (1921), charge renormalization (22), charge fluctuations (23), and nonadditive effects of dispersion interactions (24, 25). Despite these extensive studies for particles dispersed in isotropic electrolytes, the understanding of electrostatic screening in ion-doped nematic media seems to be lacking. Anisotropic electrostatic screening can potentially provide the means for controlling and engineering self-assembly of colloidal particles in nematic solvents, where a key advantage as compared with isotropic solvents may arise in defining highly anisotropic interactions and self-assembled structures. An important relevance of anisotropic electrostatic screening is also in active matter systems, as affecting major mechanisms including locomotion and energy harvesting (26). Here, we address the issue of electrostatic screening in nematic media and show that it is important by exploring an experimental model system in media with orientational elasticity.

For isotropic solvents and sufficiently dilute electrolytes, the electrostatic potential around a freely dispersed arbitrarily shaped charged particle asymptotically scales as (2729)ϕ(r)A(ψ,θ;λDI)exp (r/λDI)r,(r)(1)with r as the radial distance, ψ as the azimuthal angle, θ as the polar angle, and λDI as the constant (isotropic) Debye screening length. The anisotropy function A(ψ,θ;λDI) captures the shape effects of the particle on ϕ(r) and is independent of ψ and θ for spherical particles at any value of λDI. In general, A becomes more strongly dependent on (ψ, θ) when λDI is small, whereas A(ψ,θ;λDI)1 for λDI, meaning that ϕ(r) ∼ 1/r for any particle. In other words, anisotropies stemming from the particle shape are present even asymptotically far from the particle, contrasting the unscreened case that behaves as a point charge for r → ∞. Furthermore, note that dilute isotropic electrolytes are characterized by a single value of λDI independent of the particle orientation.

In contrast, the host material anisotropy in screened electrostatic interactions is a novel open challenge, centered at the question of how electrostatic screening is changed when the medium is characterized by a dielectric tensor rather than a dielectric constant. Spatially dependent dielectric anisotropies occur sometimes in isotropic liquids, for example, near solid-water interfaces (30, 31), but is usually confined in only a small region of space. A more general example of controllable anisotropic materials are liquid crystals (LCs), where the anisotropy is described by one dielectric coefficient along the primary dielectric tensor axis, also called the director n ≡ − n, and a second dielectric coefficient in the perpendicular direction. For a hypothetical, everywhere radial director around a colloidal sphere, only the dielectric tensor components projected in the radial direction contribute, and hence, the Bjerrum length and, consequently, the Debye screening length are renormalized with a constant, and one can just use Eq. 1 with A being constant. However, a completely different situation arises when the director configuration around the particle does not have the same symmetry as the particle itself. The simplest nontrivial example would be a constant director field along the z axis surrounding a spherical particle, and in this case, the screening will become direction dependent, as we shall see here. Furthermore, as it is the more usual case in anisotropic nematic media with dispersed particles (3236), the director field is usually spatially dependent and varies in space because of various geometries, surface effects, and external fields, leading to rich and diverse elasticity-mediated anisotropic interparticle interactions. Multipole expansions have been used to describe elasticity-mediated colloidal interactions in LCs, drawing parallels to electrostatic interactions (3739). In general, elasticity-mediated interactions in LCs are accompanied by screened electrostatic and dispersive (London–van der Waals) (40) interactions; however, the previous studies of such colloidal systems with LC hosts were done for highly anisotropic rod- and disc-shaped particles, so that the role of the anisotropy of colloidal particles and that of the LC medium were not separated because of the particle’s shape anisotropy and so far explored while probing phase behavior and self-assembly of colloidal superstructures (41, 42).

Here, we explore anisotropic colloidal interactions in electrostatically screened near-spherical charged colloids to develop a generalized understanding of electrostatic interactions in colloids, subjected to and determined by the material dielectric anisotropy. Experimentally and theoretically, we use so-called charged dumpling particles (with almost spherical shape) as a charged colloidal model system because they can become appreciably charged in a simple LC such as 5CB (4-cyano-4′-pentylbiphenyl), with a weak-enough elastic interaction that allows competition with the electrostatic forces. We calculate the effective pair interaction under the assumption that elastic, dispersive, and electrostatic interactions are additive, just like in the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory of charged-screened spheres. Furthermore, the electrostatic part is treated within linear screening theory in combination with a far-field multipole expansion approach, on the same level that typically elastic colloidal interactions are treated. Last, we compare the theoretically calculated interactions with experiments, finding good qualitative agreement.


Charged colloidal dumpling particles dispersed in a nematic electrolyte as model system

As our charged colloidal model system, we use particles (Fig. 1A) of “dumpling-like” shape. Their rough shape and overall dimensions are close to those of a sphere with a diameter 2a = 1 μm. The hydrothermal synthesis and surface treatment of these particles allow us to control the colloidal charge in a rather broad range of values without issues (typical for other particles in nematic solvents) of nonuniformity of charging. These colloidal dumplings were dispersed in 5CB at low concentration (<1000 parts per million) to obtain well-separated colloidal particles. In Fig. 1 (B to E), we show microscopy images of a single colloidal dumpling obtained in different imaging modes. The colloidal dumplings have homeotropic anchoring on their surfaces, and the symmetry of resulting director n(r) distortions around particles (Fig. 1F) is of the “quadrupolar” type, with an encircling half-integer disclination loop (“Saturn ring”) (43, 44). The in-plane diffusion of the colloidal dumplings due to Brownian motion (fig. S8, A and B, and movie S1) is anisotropic with respect to the LC far-field director n0 with diffusion coefficients D/D=1.49 to 1.54 (fig. S8C and movie S1), and this is close to theoretical predictions for spheres, D/D = 1.72 (45). We have prepared and used charged and uncharged particles in our experiments, where the effective charge was controlled in a broad range. The number of elementary charges e on the particles surface Z = 0 to 350 was determined using their electrophoretic motion between two in-plane electrodes placed perpendicular to n0 in a planar nematic cell (the uncharged particles did not move in response to applying an electric field). Dumpling charged particles were moving along n0 toward a negative electrode when a DC electric field was applied between the electrodes. The velocity of the particles depends on their charge (Fig. 1G) and the strength of the electric field. The displacement of particles was tracked using video microscopy, which allows us to estimate the effective charge Ze, from the balance of the Stokes viscous drag force and the electric force (4042). To probe only the electrostatic pair interactions between charged colloidal dumplings, we measured their pair interactions in the isotropic phase of 5CB, where the contribution due to LC elastic forces is eliminated. When brought nearby with the help of optical tweezers, colloidal dumplings repel from each other with a potential ΦE of tens of kBT. The effective charge number Z and Debye screening length can also be extracted from experimental pair interactions (Fig. 1H) using the DLVO equation (1, 2)ΦE(d)kBT=Z2λB[exp (a/λDI)1+a/λDI]2exp (d/λDI)d(2)with kBT being the thermal energy, λB the Bjerrum length, a the particle radius, λDI the (isotropic) Debye screening length, and d the center-to-center distance between particles. The effective charge numbers obtained by electrophoretic measurements (Fig. 1G) were in good agreement with values obtained from the electrostatic interaction potential (Fig. 1H). The Debye screening lengths obtained from fitting the interaction potentials were within the range of λDI = 300 to 1000 nm measured for 5CB samples in our experiments using impedance spectroscopy.

Fig. 1 Charged dumpling colloidal particles in a nematic LC as model system for anisotropic charged colloids.

(A) Scanning electron microscopy image of dumpling colloidal particles. (B to D) Polarizing and (E) bright-field microscopy textures of dumpling particles with homeotropic anchoring in a planar nematic cell between crossed (B and C) and parallel (D) polarizers P and A, with (C) and without (B) a phase retardation plate Z. (F) Schematic diagram of the director field n(r) (green lines) around a dumpling particle treated for a homeotropic anchoring: A red circle around a particle indicates a singular defect loop Saturn ring, and n0 is a far-field director. (G) Electrical effective charge of dumplings determined using electrophoresis. (H) Repulsive electrostatic pair potential ΦE between dumpling colloids in an isotropic phase of a nematic LC with λDI925 to 959 nm.

Because of the effective elastic nature of the anisotropic LC host, also uncharged colloidal particles interact via anisotropic elastic interactions (44), which, for our dumpling particles, are of quadrupolar symmetry (fig. S9) with elastic interaction potential ΦLC(d, θ) given byΦLC(d,θ)=169πKc2990 cos2θ+105 cos4θd5(3)

Note that the potential falls off as ∝1/d5 and depends on the angle θ between the uniform far-field n0 and d, which makes it strongly anisotropic, with the attraction direction at ≈40° to 50° (fig. S9, A and B). We can extract the elastic pair potential (fig. S9D) from the time-dependent separation between two attracting particles (fig. S9C), and on the basis of the elasticity measurements (using the single elastic constant K = 8 · 10−12 N), we find, for our system, the elastic quadrupole moment c = 0.1 to 0.2 μm3.

Electrostatic potential of single charged colloidal sphere in an anisotropic dielectric

Electrostatic interactions between charged particles are conditioned by the profile of the electrostatic potential surrounding the particles. This quantity exhibits anisotropic screening, as it is directly determined by the anisotropy of the host medium. We calculate this anisotropic electrostatic potential in the mean-field approach by using the Poisson-Boltzmann (PB) equation for the electrostatic potential in the nematic host with fixed director field n(r) surrounding the particle (see the Supplementary Materials)i[ϵij(r)jϕ(r)]/ϵ¯=κ2sinh [ϕ(r)],(r>a)(4)where kBTϕ(r)/e is the electrostatic potential, ϵ¯ is the rotationally averaged dielectric constant of the nematic medium, and κ1=λDI is the isotropic Debye screening length used as a “reference” decay length. Furthermore, we used the Einstein summation convention. Note that in the above, we assume that the dumpling particles can be approximated as spheres and have constant-charge boundary conditions with homogeneously distributed charge Ze and that n0 surrounding the particle has cylindrical symmetry. The host material anisotropy is given by the dielectric tensor as ϵij(r) = ϵδij + Δϵni(r)nj(r), with Δϵ = ϵ − ϵ as the dielectric anisotropy difference between dielectric tensor components projected parallel to the director ϵ and perpendicular to the director ϵ. We take, in accordance with our experiments, that ions cannot penetrate the particle with dielectric constant ϵp = 2; hence, inside the particle, one has to solve the Laplace equation2ϕ(r)=0,(r<a)(5)

The anisotropic electrostatic potential (Eqs. 2 and 3) is numerically calculated using COMSOL Multiphysics software exploiting the cylindrical symmetry, as shown in Fig. 2. We provide results for two material anisotropic regimes, one with uniform director field n = ez and second—the realistic one for our experiments—a director field with elastic quadrupole distortions n = [ez + (2czρ/r5)eρ]/∣ez + (2czρ/r5)eρ∣ (as shown in Fig. 1F). The quadrupole is derived from a multipolar expansion for the Saturn ring configuration (44).

Fig. 2 Anisotropic electrostatic potential kBT(ϕ)/re for charged spherical colloidal particle in anisotropic dielectric host.

(A and B) Potential for large (isotropic) Debye screening length and (C and D) short Debye screening length, with uniform director host c = 0 (along the z direction) and with quadrupolar director material anisotropy (c ≠ 0). (E and F) Electrostatic potential along selected directions from the particle center, at constant angle θ with respect to the z axis. Full lines correspond to uniform, and dashed lines to quadrupolar material anisotropy c = 0.2 μm3. The gray line is the isotropic electrostatic potential line. For numerical parameters, we take particle radius a = 0.5 μm with constant charge Z = 50 in a nematic electrolyte with dielectric properties ϵ = 6, ϵ = 19, and ϵ¯=10.

For the uniform director field (Fig. 2A) and a Debye screening length larger than the particle size, the diffuse screening cloud has a prolate spheroidal–like shape with the long axis coinciding with the z axis, whereas upon considering the full quadrupolar distortion, we see that the electric potential profile (i.e., the double layer) gets distorted close to the region of the Saturn ring defect (Fig. 2B) but again evolves to the prolate spheroid shape further away from the particle. Furthermore, it turns out that the spheroidal double layer has an aspect ratio for the major to minor axis equal to ϵ/ϵ. Figure 2C shows the electrostatic potential along selected directions from the particle at constant angle θ with respect to the z axis, which can be compared to the isotropic electrostatic potential, showing a similar magnitude. In Fig. 2 (D to F), the electrostatic potential for the regime of Debye screening length smaller than the particle size is shown for a uniform director (Fig. 2D) and quadrupolar director (Fig. 2E); the electrostatic potential along constant θ is shown in Fig. 2F. At these short Debye lengths, the electrostatic potential inside the particle becomes strongly inhomogeneous, and although the ions are closer to the particle, the electrostatic potential is still strongly anisotropic.

Figure 3 shows that the electrostatic potential of a spherical particle with a prolate spheroidal double layer in a nematic electrolyte behaves differently from a charged spheroidal particle with aspect ratio ϵ/ϵ in an isotropic electrolyte. Whereas the diffuse screen cloud seems to follow the shape of the spheroidal particle sufficiently close to the particle (Fig. 3, A and B), a graph on log-linear scale along a selection of cuts through the particle reveals that only the “amplitude” is anisotropic for sufficiently small λDI. However, for all bulk ion concentrations, the screening length is independent of the angle. We already alluded this well-known result in Eq. 1. In contrast, the sphere in a nematic shows—with and without topological defects—that there is an angle-dependent screening length (Fig. 3, E and F). Such an anisotropy occurs whenever the director on the particle surface (in this example, spherical) has a different symmetry than the far-field director field (cylindrical).

Fig. 3 Electrostatic potential of a charged prolate spheroid in an isotropic electrolyte compared with a charged sphere in a nematic electrolyte.

(A and B) Electrostatic potential for large and small Debye screening length, respectively, around a prolate spheroidal particle in an isotropic electrolyte with minor radius a = 0.5 μm and major radius 19/6 a. In (C) and (D), we plot the potential along selected directions from the particle center at constant angle θ with respect to the z axis on log-linear scale. In (E) and (F), we replot Fig. 2 (C and F) on log-linear scale highlighting the angle-dependent screening length found in a nematic electrolyte as opposed to an isotropic electrolyte (C and B), described by a constant Debye length. (G) Anisotropic Debye screening length λD(θ) as calculated from Eq. 10 relative to the isotropic screening length λDI in different directions.

Analytical expressions for the electrostatic potential

To determine the anisotropic electrostatic potential, we first derive an analytical expression by mapping the charged particle to an ion-penetrable spherical shell with radius R and surface charge density eσS, followed by a multipole expansion. The approach is more extensively explained in Materials and Methods, the Supplementary Materials, and (46).

The calculated electric potential isϕ(r)=αγ2ZλBIϵ2ϵ/ϵ¯3[Gm(r)+(γa)26Gq(r)+(γa)4120Gh(r)+](6)as expressed with multipolar basis functions Gi(r) (i = m, q, h, . . ) (see the Supplementary Materials for their expressions). The higher-order multipoles become more important at higher salt concentrations (smaller λDI), and the parameters α and γ also depend on a/λDI. The salt-dependent parameters α = σS/σ > 1 and γ = R/a < 1 can be determined using a fit to the numerically obtained surface potential. By fitting only the surface potential, it is found that the integral expression that is obtained before performing the multipole expansion, is practically numerically exact for r > a in a wide range of salt concentrations (see the Supplementary Materials).

As is usual for multipole expansions, Eq. 6 fails at short distances, but at large-enough r = ∣r∣, it captures the proper angle dependence, given that enough multipoles are taken into account (see the Supplementary Materials). For example, up until hexadecapolar order, Eq. 6 is numerically exact up until a/λDI2 (<5% deviation), while at higher salt concentrations, truncation at the hexadecapolar order turns out to be not sufficient. As an example, for a/λDI5, we see, even asymptotically far from the particle, that the deviation is 10 to 60%, depending on the angle with the director (see the Supplementary Materials). Last, Z is the actual charge number for ∣ϕ(r)∣ ≪ 1, but for high electrostatic potentials, Z should be interpreted as a renormalized charge density, similar to what is known in “isotropic” charged colloids (22, 47).

From Eq. 6, we derive the asymptotic scalingϕ(r)ZA(θ)exp [r/λD(θ)]r/λB(θ),(r)(7)with the anisotropy function for uniform director fieldsA(θ)=αγ2{1+(γκa)26ϵ¯3(cos2θϵ+sin2θϵ)(cos2θϵ2+sin2θϵ2)+O[(κa)4]}(8)angle-dependent Bjerrum lengthλB(θ)λBI=ϵ¯ϵ(ϵΔϵcos2θ)(9)and angle-dependent Debye screening lengthλD(θ)λDI=ϵϵϵ¯(ϵΔϵcos2θ)(10)with, for Δϵ > 0 (Δϵ < 0), a maximum (minimum) at ϵ/ϵ¯ and a minimum (maximum) at ϵ/ϵ¯. Note that λD(θ)/λDIλBI/λB(θ), as one would maybe naively think based on λDI=(8πλBIρs)1/2. We plotted Eq. 10 in Fig. 2G for Δϵ > 0, and from this, the shape of the double layer can be understood. Had we taken Δϵ < 0, the double layer would have had an oblate spheroidal shape.

Equation 7 reveals why electrostatic screening in an anisotropic medium is markedly different from screening of anisotropic particles in an isotropic medium. For r → ∞, the particle anisotropy can be fully accounted for with the anisotropy function A(θ) (see Eq. 1), whereas anisotropic media are described not only by an anisotropy function (even for spheres with nonvanishing volume) but also by an angle-dependent Debye and Bjerrum length (see Eqs. 9 and 10, respectively). The anisotropic screening length is highlighted in Fig. 3 (E and F) and is the relevant decay length even in the presence of defects. Furthermore, Fig. 3 (E and F) shows, in accordance with Eq. 7, that the anisotropy function is only relevant for a sufficiently large κa, giving rise to larger amplitude differences of the electrostatic potential for varying θ and at fixed r. Furthermore, Fig. 3F suggests that defects alter the form of A(θ), but not that of λD(θ), which is determined only by the symmetry of the far-field director. Therefore, we hypothesize that defects asymptotically behave as if they change the particle anisotropy but not the medium anisotropy.

Last, while assuming that the strength and nature of surface boundary conditions on particle surfaces do not change with adding ions, we evaluate the salt-dependent renormalization parameters α and γ for some values of a/λDI (more information in the Supplementary Materials). For a/λDI=0.5, we find γ = 1 and α = 1.06; for a/λDI=1, we find γ = 0.99 and α = 1.2; and for a/λDI=2, we find γ = 0.988 and α = 1.43. With these quantities, we can use Eq. 6 in the relevant experimental parameter regime, and the same parameters enter the expression for the effective far-field pair potential (see below).

Anisotropic electrostatic and elastic pair interactions

Charged colloidal dumplings interact in the nematic LC (Fig. 4A and movie S2), both via anisotropic electrostatic and elastic interactions, each with a different anisotropic profile. Generally, the electrostatic interaction is repulsive, whereas the nematic elastic interaction has regions (directions) of attraction and regions of repulsion. If the charge of the particles is high enough (∼300e), then the electrostatic repulsion can counterbalance the elastic attraction at any angle and the two colloidal particles stay separated from each other without getting into the full contact (Fig. 4) at any θ. The orientation of the separation vector between particles is primarily ≈55° (Fig. 4B) relative to n0. We note that for this histogram, we mapped the angles to the first quadrant because of the quadrupolar symmetry of the system. Histograms of center-to-center distance between particles (Fig. 4C) show that, for fixed angles θ ≈ 0 and 90 where we used laser tweezers to release the particle pair along these angles, the separation is d ≈ 6 and 5 μm, respectively, which is consistent with the anisotropy of the calculated electrostatic repulsion force and overall colloidal interactions in this system (Fig. 4E). This difference in the steady-state d depends on the position of the two particles relative to n0, which can result from the anisotropy in the charge distribution around colloidal dumplings (Fig. 2). Histograms of the angle θ and an overlaid polar plot of separations between two freely moving particles (Fig. 4, B and D) show that there is a preferred orientation for a separation d at ≈55° (again mapping all the data to the first quadrant), which indicate an equilibrium distance and orientation of a particle pair resulting from a competition of anisotropic elastic attraction and anisotropic electrostatic repulsion (Fig. 4, D to G), as we shall explore theoretically below. On the other hand, if the electric charge at the dumpling surface is sufficiently small (∼150e), then the elastic attraction is dominant and two colloidal particles attract and get to the full surface-to-surface contact as elastic quadrupoles (fig. S10), similar to colloidal dumplings without charge (fig. S9 and movie S3).

Fig. 4 Pair interactions of charged dumpling particles in a nematic LC.

(A) Bright-field micrograph of two interacting colloidal dumplings. (B) Histogram showing the preferred angle θ between the center-to-center separation vector and the far-field director for a pair of two highly charged particles. (C) Histograms showing preferred separations between two highly charged particles (∼300e per particle) when their separation vector is at θ ≈ 0 (blue) and θ ≈ 90 (red). (D) Analytically calculated total effective potential as the sum of contributions from (E) anisotropic electrostatic interactions, (F) anisotropic elastic interactions, and (G) van der Waals dispersion interactions. The interaction potentials are calculated for particle radius a = 0.5 μm, particle charge number Z = 750, elastic quadrupole moment c = 0.1 μm3, Hamaker constant AH = 1 kBT, and nematic electrolyte with isotropic Debye screening length λDI=0.5μm. Black crosses in (D) show the overlaid angular dependent separations between two interacting charged particles: Particles do not come in direct contact at any angle. Note that experimental preferred angles and separations are consistent with the calculated competing electrostatic and elastic potentials, with the minimum of the total effective interaction potential at ∼55°.

The effective pair interaction between anisotropic charged colloids is determined theoretically by splitting—in the spirit of the DLVO theory—the total interaction potential Φ(d, θ), as the sum of screened electrostatic ΦE(d, θ), van der Waals ΦvdW(d), and, because we are in a nematic host, the nematic elastic interactions ΦLC(d, θ)Φ(d,θ)=ΦE(d,θ)+ΦvdW(d)+ΦLC(d,θ)(11)

The van der Waals interaction can be derived using Hamaker–de Boer theoryΦvdW(d)=AH3[a2d24a2+a2d2+12ln (14a2d2)](12)which fails at center-to-center distances close to contact (d ≈ 2a), where the (quantum-mechanical) Born repulsion becomes important, and for large d, where relativistic effects become important. For this interaction, the anisotropy enters only the Hamaker constant, but not in the coordinate-dependent part of the expression, assuming the nonrelativistic limit (48). Last, the elastic interaction ΦLC(d, θ) is given by the quadrupolar far-field expression Eq. 3.

We determine the effective anisotropic screened electrostatic interaction ΦE(d, θ) from the asymptotic expression of the electrostatic potential (Eq. 6) within the linear superposition approximation (LSA; for full derivation, see the Supplementary Materials). The anisotropic screened electrostatic interaction readsΦE(d,θ)kBT=α2γ4Z2λBIϵ2ϵ/ϵ¯3[Gm(d,θ)+(γa)23Gq(d,θ)+245(γa)4Gh(d,θ)+](13)now treated on the same level of (multipolar) approximation as the elastic interactions, where Gm, q, h(d, θ), α, and γ are the same as the ones derived for the single-particle case (see Eq. 6). Note also that in the limit of low salt concentration, only the monopole term is relevant, and the above equation reduces to the anisotropic Yukawa formΦE(d,θ)kBTα2γ4Z2λB(θ)exp [d/λD(θ)]d(14)which is reminiscent of the well-known standard isotropic DLVO potential (Eq. 2).

In Fig. 4D, we show the total analytically calculated interaction potential for parameters, in qualitative range of our experiments. The position of a local minimum of several kBTs is found at θ ∼ 55, which is in good qualitative agreement with experiments (Fig. 4B). The position of the theoretically calculated local minimum is also in good agreement with the overlaid experimentally determined orientations of the separation vector (Fig. 4D), where a pair of interacting particles reside primarily also at ≈55°. In Fig. 4E, we show the contributions of the anisotropic electrostatic effective interaction, which also shows a good agreement with the experiment (Fig. 4C); in Fig. 4F, the elastic interaction; and in Fig. 4G, the van der Waals interaction. The van der Waals interaction is of shortest range and does not contribute substantially to the total interaction potential, whereas elastic and screened electrostatic potentials clearly compete, and it is their detailed balance that determines the overall interparticle potential. Last, note that if we had used the isotropic DLVO potential (Eq. 2) to calculate the total potential (and using same material parameters), the predicted local minimum would shift to an angle θ ∼ 49, which underlines the clear role of the electrostatic anisotropy.

The total interaction potential (Eq. 11) exhibits a range of possible qualitatively different interparticle interaction regimes, as we show in Fig. 5, which depend on the particle charge Z and the host electrolyte screening length λDI (salt concentration). Notably, we calculate the total regimes for a smaller elastic interaction c = 0.02 μm3 than in our experiments (c ∼ 0.2 μm3), which makes the anisotropic electrostatic DLVO-type interaction of more similar magnitude as the nematic elasticity at the reported particle charges, which are lower than the one chosen in Fig. 4D. We want to stress, however, that the electrostatics is based on a far-field multipole expansion in combination with the LSA; hence, our theory underestimates the repulsion when the double layers of the particles overlap at sufficiently low salinity and low particle separations (fig. S5), as is also common in isotropic DLVO theory (5), not to mention because of the currently unknown nonadditive effects between elasticity and electrostatics, or close-approach elastic effects.

Fig. 5 Effect of particle charge and salt concentration on total interaction pair-potential of spherical charged colloidal particles with quadrupolar nematic dielectric anisotropy.

(A to I) The total potential is calculated for particle radius a = 0.5 μm and elastic quadrupole magnitude c = 0.02 μm3 for varying particle charge Z and different screening lengths λDI. Note the difference in color scales for different screening lengths. Black dots in (E) show experimental results, i.e., relative particle separations, overlying the theoretically calculated interaction potential in the relevant material parameter regime (Z = 150, λDI=0.5μm). Note that particles do not come in direct contact at any angle.

In Fig. 5 (A to C), we show the situation when the screening length is larger than the particle size. For low-enough charges, a global minimum located at approximately 49 is separated from a shallow local minimum with depth ∼0.001kBT. Both local and global minimum disappear when the particle charge is increased (Fig. 5, B and C), and the interaction potential is repulsive for all distances and directions. For Debye lengths similar to the particle size, we see that higher charges are needed to overcome the attractive elastic interaction. At the same particle charge but smaller screening length (Fig. 5D), the electrostatic interaction is too weak to overcome the elastic interaction, resulting in a purely attractive direction in the interaction potential. Such a situation is reminiscent of what we also observed experimentally (but for slightly different exact parameters) in fig. S10, where a low charge results in particle coagulation. When increasing Z further, again a local minimum (Fig. 5E) that is deeper than the low-screening case emerges. This local minimum is in good agreement with experimentally determined angular dependent separations between interacting particles overlaid over the calculated data, even though we choose a smaller value for the elastic quadrupole than was experimentally measured (fig. S9); therefore, a lower particle charge is needed to generate the metastable local minimum. The depth of this minimum becomes smaller when Z increases even further (Fig. 5F). The trend that deeper attractive minima can be attained for larger salt concentrations, but that higher charges are needed, is something that we also see for double layers smaller than the particle size (Fig. 5, G to I). When comparing Fig. 5 with Fig. 4, we see that increasing the strength of elastic interaction at the expense of higher Z also gives rise to deeper attractive minima.


Summarizing, we introduced a screened electrostatic colloidal model system that can become appreciably charged in a nematic LC, with particles of almost spherical shape. Theoretically, we discussed that dielectric anisotropy of the colloidal host, given in nematic fluids by the director field, gives rise to anisotropic screening of the electrostatic potential when there is a mismatch between the nematic and particle symmetry, and, in turn, also to the electrostatic effective pair interaction. In our system of nematic colloids, the screened electrostatic interaction is inherently combined with effective nematic elastic interactions, which leads to different interaction regimes where particles are (i) repelled for all distances and angles, (ii) are subjected to a weak local minimum of <1kBT such that they can still move because of thermal fluctuations (Fig. 4, B and D), and (iii) are dominated by the elastic interactions with distinct attractive and repulsive directions (fig. S10). In our experiments with charged colloidal dumplings, we have observed regime (ii) at Ze ∼ 300e and regime (iii) at Ze ∼ 150e, whereas to access for regime (i), we would need to achieve even higher particle charges and/or more deionized samples. However, we note that regime (i) was achieved in experiments with highly shape-anisotropic particles in (41), where the small diameter of rod-like colloidal particles allowed for fully overpowering elastic interactions by the electrostatic ones. Our model is consistent with the notion that controlling particle dimensions provides the means of shifting the balance between the elastic and electrostatic forces at given colloidal surface charge and ionic conditions.

We have shown, within an experimentally accessible parameter range, that even effectively spherical particles exhibit anisotropic electrostatic interactions in LCs. Experimental results indicate that both elastic and electrostatic interactions of charged particles in nematic LCs are relatively long ranged and anisotropic with respect to the director, so that the colloidal behavior depends on their interplay. While charging could be controlled from neutral to hundreds of elementary effective charges per single particle, we could show that the interparticle forces could be dominated by elasticity or by electrostatics in the two limiting regimes, with both elastic and electrostatic interactions being highly anisotropic. While we focused on thermotropic nematics with accessible range of host medium’s Debye screening lengths in the range of 300 to 1000 nm and on colloidal particles with the accessible range of surface charges (0 to 350)e, our findings do indicate a plethora of colloidal behaviors arising from the interplay of electrostatic and elastic interactions with salient anisotropic features, consistent with theoretical modeling.

This study can be, in the future, extended further by exploiting the elasticity and electrostatics interplay for lyotropic water-based LC colloidal systems, where Debye screening lengths can be much shorter, on the order of several nanometers, as well as highly deionized LCs that potentially could allow for accessing the range of Debye screening lengths from several nanometers to 10 μm. As shown in Fig. 5, tuning the screening length can change the relative position of a local minimum in the effective pair potential. Moreover, we note that lyotropic systems might have additional features that our theory does not explore yet, being often a five-component mixture (an isotropic solvent such as water, LC particles, cations, anions, and the larger colloidal particles), compared with the thermotropic four-component suspension in this paper (LC, cations, anions, and colloidal particles). Note that lyotropic systems can also have a dielectric anisotropy that couples the director with electrostatics (49), just as in the thermotropic systems under consideration here. Furthermore, when the LC lyotropic particles are smaller colloidal (nano)particles, we envisage tuning of the dielectric, elastic, and possibly flexoelectric properties of the nematic host medium by changing the particle functionality. Tuning of elastic properties by charged nanoparticles as a function of particle charge and salt concentration has already been explored theoretically (50). Moreover, lyotropic systems can be made active, giving rise to possibly new unexplored hydrodynamic-electrokinetic active processes, which may be interesting also in a biological setting.

From the point of view of particle designs, these studies could be extended to patchy particles with different densities or even signs of charges, potentially allowing for different electrostatic multipoles, whereas our study showed that homogeneously charged spheres only give rise to even anisotropic Yukawa multipoles. Furthermore, while highly anisotropic disc- and rod-shaped particles have already been studied (4042), there is a considerable range of possibilities in defining colloidal behavior also by the particle shape and surface treatment for different boundary conditions for the LC director at colloidal surfaces.

As further main theoretical results, we derived asymptotic theoretical expressions for the electrostatic potential and the resulting pair interaction for homogeneous director configurations, which highlights an anisotropy not only in the screening length but also in the prefactor, where the latter also occurs for anisotropic particles in isotropic solvents. Both anisotropies together predict that local minima in the total pair potential are shifted (Fig. 4, B and D) in terms of the equilibrium angle compared with an isotropic electrostatic interaction, in line with our experimental observations. Last, one obvious extension to the theory is to numerically investigate the effect of close distances between particles and relaxing the requirement of additivity in the pair potential by coupling the full Landau–de Gennes theory in terms of the tensor order parameter with electrostatics. We will explore this in future work.

More generally, our work contributes to the generalization and extension of the DLVO interactions to ubiquitous anisotropic soft matter systems, such as complex fluids and anisotropic colloids. While nematic colloids formed by a thermotropic LC host and near-spherical colloidal inclusions provide validation of our theoretical findings, the concepts introduced here can be applied in biological contexts of highly structured biological cell interior and membranes, active matter systems with the additional forms of anisotropy stemming from activity, lyotropic LCs with varying degrees of orientational and partial positional ordering, ionic fluids, and so on. While the experiments and model focused on even anisotropic Yukawa multipoles formed by like-charged spherical particles, future studies can extend our concepts to odd anisotropic Yukawa multipoles via a heterogeneous surface charge distribution or nonspherical particle shape on the electrostatic side of the spectrum, and to elastic monopoles through hexadecapoles and higher-order multipoles on the elastic side. It will be of interest to consider further the (nonadditive) effects of various topological defects on counterion distributions, the role of flexoelectricity and surface polarization, surface anchoring effects, charge regulation, flow (51), and how similar concepts work in LC mesophases with different point group symmetries and partial positional ordering. Specifically, one can expect to have the double layer–driven realignment of the nematic and even electrostatically controlled surface anchoring transitions as function of salt concentration and interparticle distance in the regime where the dielectric torque is of similar magnitude as the elastic torque (42, 5254). Furthermore, flexoelectric topological defects can substantially alter charge distributions on the colloidal particles (55). These are just a few of the many examples of additional higher-order phenomena of which the effects on the effective interparticle potential are left for future studies. Overall, our findings will contribute to the soft matter toolkit for forming colloidal composite materials with pre-engineered structure and composition of the constituent building blocks and will also help further in understanding nonequilibrium phenomena in these systems (56).


Synthesis and characterization of charged colloidal dumplings

The dumpling colloidal particles (Fig. 1A) were prepared using the hydrothermal synthesis method as reported earlier (57). The chemical ingredients used for synthesis, ytterbium nitrate hexahydrate [Yb(NO3)3 6H2O], yttrium nitrate hexahydrate [Y(NO3)3 6H2O], erbium nitrate pentahydrate [Er(NO3)3 5H2O], and sodium fluoride (NaF) were all purchased from Sigma-Aldrich. Sodium hydroxide (NaOH) was purchased from Alfa Aesar. Octanoic acid (OA) was purchased from Acros Organics. In a typical synthesis, 130 mg of Y(NO3)3, 40 mg of Yb(NO3)3, and 10 mg of Er(NO3)3 were mixed with 10 ml of deionized water and 13.5 ml of ethanol. After forming a clear transparent solution, 0.35 g of NaOH and 1.83 g of OA were added into the above solution and stirred continuously at 50C for 30 min. Then, 9 ml of 0.2 M NaF solution in deionized water was added dropwise to the above solution and stirred continuously for 30 min at 50C. The mixture was transferred to a 40 ml Teflon-lined autoclave and kept in an oven at 200C for 7 hours. After the reaction, the autoclave was allowed to cool down to room temperature naturally. The particles precipitated at the bottom of the reaction vessel were collected by centrifugation, washed with ethanol and deionized water in sequence, and, lastly, dispersed in 5 ml of cyclohexane. The reaction yields OA functionalized dumpling-shaped particles with an average size of 1 μm, as demonstrated by the scanning electron micrograph of the particles deposited onto a silicon substrate (Fig. 1A). To induce positive surface charges, the particles were treated with an acidic solution. Briefly, 200 μl of concentrated hydrochloric acid (HCl) was mixed with 5 ml of deionized water, and 2.5 ml of the particle solution in cyclohexane was added dropwise to the acid solution and stirred continuously for 12 hours. During this reaction, the OA molecules originally attached to the particle’s surface got detached, leaving the particle positively charged. The uncapped particles were collected by the centrifugation and subsequently dispersed in 2.5 ml of ethanol for further use. To prepare the colloidal dispersions, particle solution was mixed with 5CB (Frinton Laboratories), followed by solvent evaporation at 70C for 2 hours and cooling to nematic phase under vigorous mechanical agitation. The nematic dispersions were infiltrated to a glass cell by capillary action and sealed using a fast-setting epoxy. We used planar nematic cells with tangential boundary conditions for experimental observations and video tracking of anisotropic pair interactions of charged and uncharged particles in nematics. To promote unidirectional tangential boundary conditions at the substrate interface, the inner surfaces of the glass substrates were coated with polyvinyl alcohol and then rubbed unidirectionally.

We used an experimental setup built around an inverted Olympus IX81 microscope and a 100× (numerical aperture 1.4) oil objective to perform bright-field and polarizing microscopy observations. Translational and rotational motion of colloidal particles was recorded with a charge-coupled device camera (Flea, Point Grey) at a rate of 15 frames per second, and the exact x-y position (fig. S8 and movie S1) of the dumpling particles as a function of time was then determined from captured sequences of images using motion tracking plugins of ImageJ software. Optical manipulation of dumpling particles was realized with a holographic optical trapping system operating at a wavelength of 1064 nm and integrated with our optical microscopy system. To measure the Debye screening length of the 5CB samples, we used impedance spectroscopy with the Schlumberger SI 1260 impedance analyzer.

Calculation of electrostatic potential

To calculate the electrostatic potential in an anisotropic dielectric medium, we start from Gauss’ law, which is given inside a spherical particle with radius a due to the absence of free charges by·D(r)=0,(r<a)(15)with D(r) = ϵ0ϵpE(r) as the displacement field expressed in terms of the vacuum permittivity ϵ0, particle dielectric constant ϵp, and electricfield E(r). Outside the particle, in the nematic, we have ions with number densities ρ±(r), and hence, the Gauss law reads (in SI units)·D(r)=e[ρ+(r)ρ(r)],(r>a)(16)where now, Di(r) = ϵ0ϵijEj(r), with ϵij(r) as the components of the (symmetric) dielectric tensor. We can express the electric field in terms of the electrostatic potential kBTϕ(r)/e, and using the result that within the mean-field approximation the ion densities are Boltzmann distributed, we find Eqs. 4 and 5 of the main text, solved numerically with COMSOL Multiphysics under the assumption of a constant charge density eσ. In the Supplementary Materials, we derive Eqs. 4 and 5 also from a free-energy approach.

For ∣ϕ(r)∣ ≪ 1 and a constant dielectric tensor, Eq. 4 simplifies to2ϕ(r)=0,(r<a)(17)and[(ϵij/ϵ¯)ijκ2]ϕ(r)=0,(r>a)(18)to be matched by the constant-charge boundary conditionν̂i[ϵijjϕ(r)r=a+ϵpiϕ(r)r=a]/ϵ¯=4πλBIσ(19)with ν̂i being an outward-pointing unit normal vector. To obtain analytical solutions is, however, difficult because of the different symmetry inside the particle compared with outside the particle, which prevents us to find a closed-form expression for ϕ(r) while satisfying the constant-charge Neumann boundary condition. We can, however, find a very accurate analytical approximation used in the main text, for which we will sketch the approach here and leave the details of the calculations for the Supplementary Materials.

The most important step in finding an analytical solution is to observe that an approximate solution can be found by solving the auxiliary problem of an ion-penetrable charged shell with surface charge density eσS and radius R[(ϵij/ϵ¯)ijκ2]φ(r)=4πλBIσSδ(rR)(20)where σS and R need to be fitted to the numerically obtained electrostatic potential. Then, it turns out thatϕ(r)raφ(r)ra(21)

In general, and especially at high salt concentrations, Ra and σS ≠ σ, except in the limit where κ → 0. The advantage of the auxiliary problem is that the solution has a closed-form expression with only a double integral left to be evaluated, in terms of the analytically known anisotropic Debye-Hückel (DH) Green’s function G(r, r′), defined by[(ϵij/ϵ¯)ijκ2]G(rr)=4πδ(rr)(22)

Hence, we only have to determine the parameters γ = R/a and α = σS/σ based on a two-parameter fit of the numerically obtained surface potential to get the electrostatic potential everywhere outside of the particle. The resulting integral expression of φ(r) turns out to be almost indistinguishable from the real ϕ(r) (see the Supplementary Materials for comparative figures). Evaluating a double integral is computationally less expensive than solving the (linearized) PB equation, but the real value of the integral representation comes when calculating pair interactions (see next subsection).

Moreover, the integral expression gives analytical insight. It is now possible to perform a multipole expansion because we have an integral representation of the electrostatic potential outside the particle, as a convolution of a singular charge distribution with the anisotropic DH Green’s function. Performing this expansion, the remaining double integrals can be evaluated to find that the decay length is given by Eq. 10. As is common with multipole expansions, they fail at short distances from the particle, as can be seen from the comparative figures supplied in the Supplementary Materials, but still capture the correct angle dependence for sufficiently large distances.

Calculation of the screened electrostatic pair interaction potential

Within linear screening theory ∣ϕ(r)∣ ≪ 1, ion entropy terms do not contribute to the effective pair potential, and hence, the electrostatic part of the pair potential is given byΦEkBT=12drq(r)ϕ(r)2UselfkBT(23)with q(r) = Σi = 1,2σδ(∣rRi∣− a), with R1 and R2 as the center-of-mass position of particles 1 and 2, respectively, and Uself is the self energy of a single particle. Now, ϕ(r) is the dimensionless electrostatic potential of the two-body problem. Using the LSA, which entails that the two-body electrostatic potential is given by the sum of the single-particle electrostatic potentials of the two particles, gives the DLVO expression (Eq. 2) if one uses the stress-tensor method (46). However, applying the LSA directly to Eq. 23 gives the wrong result because it inappropriately accounts for ion exclusion from the hard core of one particle caused by the double layer of the other particle, which is a curious peculiarity of the free-energy route to pair interactions in the theory of charged colloids (28).

Using the spherical shell renormalization method, on the other hand, we find that the effective pair potential can be approximated asΦEkBTdr qS(r)φ2S(r)2USselfkBT(24)

Here, qS(r) = Σi = 1,2σSδ(∣ rRi∣− R) is the charge distribution of two ion-impenetrable charged shells, with the same center-of-mass positions as the spherical particles. Note that singular-charge distributions have an infinite self-energy USself that we have to subtract in Eq. 24 by using an appropriate regularization procedure. For example, one way is by giving the shells a finite thickness and taking the thickness in the final step of the calculation to zero. It can straightforwardly be shown that the two-shell electrostatic potential is simply φ2S(r) = Σi = 1,2φ(rRi). Therefore, the real benefit of this method is that Eq. 24 is determined by the same α and γ that are determined from the single-particle problem. This method resembles how the electrostatic part of the DLVO expression can be obtained by solving the auxiliary problem of a spherical shell or a point charge, accounting properly for ion-hard core exclusions, although still using the free energy route. The point/shell charge value together with the ions contained within r < a equals the charge on the particle [see, for a more detailed discussion, the appendix in (28) for the shell calculation and (46) for the point-charge method]. Unfortunately, a similar method to determine σS does not apply here since ϕ(r) is inhomogeneous for r < a when the particle is dispersed in an anisotropic dielectric medium, which therefore requires a numerical fit. However, using the shell method does give the correct expression within LSA via a free-energy route.

Furthermore, the DLVO theory is just the first-order term in a complicated series expansion, where the higher-order terms become more important in the case of high salt concentrations and strong double-layer overlap beyond the LSA (2). We can therefore expect that Eq. 24 becomes progressively more inaccurate at high salt concentrations and strong double layer overlaps (5). Comparing with numerical calculations, we show in the Supplementary Materials that this is indeed the case.

To evaluate Eq. 24, one has to perform four integrals numerically, which is computationally less expensive than solving the (three-dimensional) PB equation, followed by a stress tensor or free-energy integration. To progress further, one can perform again a multipole expansion, but this time of the pair potential from the shell method to obtain an analytically tractable expression that gives more insight in the physics of the effective pair interaction in anisotropic media. See the main text Eq. 13 and the Supplementary Materials for the derivation of the multipole expansion, as well as comparison with numerical calculations of the full nonlinear theory (but at low charges). Second, we choose to use the multipolar expansion in the main text because we want to treat the screened electrostatic interaction on the same footing as the elastic interaction, which is given on the level of a multipolar expansion as well.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: J.C.E. acknowledges fruitful discussions with S. Čopar and A. Šarlah. Funding: J.C.E. acknowledges financial support from the European Union’s Horizon 2020 programme under the Marie Skłodowska-Curie grant agreement no. 795377 and from the Polish National Agency for Academic Exchange (NAWA) under the Ulam Programme grant no. PPN/ULM/2019/1/00257. M.R. acknowledges financial support from the Slovenian Research Agency ARRS under contracts P1-0099, J1-1697, and L1-8135. B.S., H.M., and I.I.S. acknowledge funding from the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under award ER46921, contract DE-SC0019293 with the University of Colorado Boulder. Last, the authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program (The Mathematical Design of New Materials) when work on this paper was undertaken. I.I.S. also acknowledges the hospitality of the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara, where he was working on this publication during his extended stay and where his research was supported, in part, by the U.S. NSF under grant no. NSF PHY-1748958. This work was supported by EPSRC grant number EP/R014604/1. Author contributions: J.C.E. performed the numerical and analytical theoretical calculations under the supervision of M.R. B.S. and H.M. performed experiments and analyzed experimental data under the supervision of I.I.S. All authors contributed to writing and discussing the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article