## Abstract

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.

## INTRODUCTION

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 (*4*–*7*), ion transport (*8*), the influence of external fields (*9*–*11*), and predicting phase behavior, such as demixing (*12*, *13*) and crystallization (*14*, *15*). Furthermore, there are many theoretical studies on many-body effects (*16*–*18*), charge regulation (*19*–*21*), 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 (*27*–*29*)*r* as the radial distance, ψ as the azimuthal angle, θ as the polar angle, and **r**) and is independent of ψ and θ for spherical particles at any value of **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

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 (*32*–*36*), 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 (*37*–*39*). 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.

## RESULTS

### 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 2*a* = 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 **n**_{0} 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 **n**_{0} in a planar nematic cell (the uncharged particles did not move in response to applying an electric field). Dumpling charged particles were moving along **n**_{0} 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 (*40*–*42*). 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 *k*_{B}*T*. 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*)*k*_{B}*T* being the thermal energy, λ_{B} the Bjerrum length, *a* the particle radius, *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

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

Note that the potential falls off as ∝1/*d*^{5} and depends on the angle θ between the uniform far-field **n**_{0} 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 μm^{3}.

### 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)*k*_{B}*T*ϕ(**r**)/*e* is the electrostatic potential, *Ze* and that **n**_{0} surrounding the particle has cylindrical symmetry. The host material anisotropy is given by the dielectric tensor as ϵ* _{ij}*(

**r**) = ϵ

_{⊥}δ

*+ Δϵ*

_{ij}*n*(

_{i}**r**)

*n*(

_{j}**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 equation

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** = **e*** _{z}* and second—the realistic one for our experiments—a director field with elastic quadrupole distortions

**n**= [

**e**

*+ (2*

_{z}*cz*ρ/

*r*

^{5})

**e**

_{ρ}]/∣

**e**

*+ (2*

_{z}*cz*ρ/

*r*

^{5})

**e**

_{ρ}∣ (as shown in Fig. 1F). The quadrupole is derived from a multipolar expansion for the Saturn ring configuration (

*44*).

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

### 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*G _{i}*(

**r**) (

*i*= m, q, h, . . ) (see the Supplementary Materials for their expressions). The higher-order multipoles become more important at higher salt concentrations (smaller

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

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

### 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 (∼300*e*), 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 **n**_{0}. 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 **n**_{0}, 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 (∼150*e*), 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).

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*, θ)

The van der Waals interaction can be derived using Hamaker–de Boer theory*d* ≈ 2*a*), 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*G*_{m, 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

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 *k*_{B}*T*s 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 *c* = 0.02 μm^{3} than in our experiments (*c* ∼ 0.2 μm^{3}), 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.

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.001*k*_{B}*T*. 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.

## DISCUSSION

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 <1*k*_{B}*T* 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* ∼ 300*e* and regime (iii) at *Ze* ∼ 150*e*, 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 (*40*–*42*), 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*, *52*–*54*). 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*).

## MATERIALS AND METHODS

### 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(NO_{3})_{3} 6H_{2}O], yttrium nitrate hexahydrate [Y(NO_{3})_{3} 6H_{2}O], erbium nitrate pentahydrate [Er(NO_{3})_{3} 5H_{2}O], 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(NO_{3})_{3}, 40 mg of Yb(NO_{3})_{3}, and 10 mg of Er(NO_{3})_{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 50^{∘}C 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 50^{∘}C. The mixture was transferred to a 40 ml Teflon-lined autoclave and kept in an oven at 200^{∘}C 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 70^{∘}C 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}ϵ_{p}**E**(**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 _{i}*(

**r**) = ϵ

_{0}ϵ

*(*

_{ij}E_{j}**r**), with ϵ

*(*

_{ij}**r**) as the components of the (symmetric) dielectric tensor. We can express the electric field in terms of the electrostatic potential

*k*

_{B}

*T*ϕ(

**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 to**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*_{S} and *R* need to be fitted to the numerically obtained electrostatic potential. Then, it turns out that

In general, and especially at high salt concentrations, *R* ≠ *a* 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

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*q*(**r**) = Σ_{i = 1,2}σδ(∣**r** − **R*** _{i}*∣−

*a*), with

**R**

_{1}and

**R**

_{2}as the center-of-mass position of particles 1 and 2, respectively, and

*U*

_{self}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

Here, *q*_{S}(**r**) = Σ_{i = 1,2}σ_{S}δ(∣ **r** − **R*** _{i}*∣−

*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

_{2S}(

**r**) = Σ

_{i = 1,2}φ(

**r**−

**R**

*). 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*

_{i}*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 MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/5/eabd0662/DC1

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.

## REFERENCES AND NOTES

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

- Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).