Electrostatic interactions between spheroidal dielectric particles

Theory is developed to address the signiﬁcant problem of electrostatic interactions between charged polarisable dielectric spheroids. The electrostatic force is deﬁned by particle dimen-sions and charge, dielectric constants of the interacting particles and medium, and interparticle separation distance; and it is expressed in the form of an integral over the particle surface. The switching behaviour between like charge repulsion and attraction is demon-strated as depending on the ratio of the major and minor axes of spheroids. When the major and minor axes are equal, the theory yields a solution equivalent to that obtained for spherical particles. Limiting cases are presented for non-polarisable spheroids, which describe electrostatic behaviour of charged rods, discs and point charges. The developed theory represents an important step towards comprehensive understanding of direct interactions and mechanisms of electrostatically driven self-assembly processes.


I. INTRODUCTION
Direct interactions and electrostatic forces often serve as a basis for novel self-assembly mechanisms, where the interacting particles combine to form larger ordered structures, typically when subjected to an external stimulus (solvent polarity, pH factor, irradiation, temperature) and driven by thermodynamic and other constraints. Significant advances 1 have been reported on designing nanoparticles with specific shapes, morphological features and interfaces that result in directional interactions in order to achieve the desired extended structures and their functionalities. Breakthroughs in particle synthesis led to the production of particles in the shape of rods, 2 cones 3 and discs, typically containing silica, metals, metal oxides, [4][5][6] and polymers, 7 with high yield and size/shape selectivity; these include some elegant examples of rods and ellipsoids of Au-Pt 8 , CdSe 9 , gold 10 , gibbsite 6 , and polymer latex. 11 These new approaches to particle synthesis have offered a diverse spectrum of particle anisotropy and clustering behaviour, including the formation of low symmetry clusters 12 and spherical self-assembled objects, 13 chain-like structures 13 and bundling. 14 Equilateral polygonal platelets have been lithographically fabricated to demonstrate that colloidal interactions and self-assembly in anisotropic nematic fluids can be effectively tailored through the control over the particles' shapes. 15 Some additional chemical and biological application areas reliant on the accurate description of electrostatic interactions between objects with spheroidal, or near spheroidal, shapes are fullerenes of higher order (e.g. C 70 ), 16 complex polyoxometalates (POMs) a) Author to whom correspondence should be addressed. Electronic mail: Elena.Besley@nottingham.ac.uk (e.g. Preyssler-type POMs), 17 elliptocytes (abnormally shaped red blood cells), 18 and some proteins. 19,20 Moreover, non-sphericity affects the self-assembly of many other types of nanoparticles, 21 the formation photonic and liquid crystals, 22,23 and light scattering. 24 Therefore, it is crucial to understand the correlation between the shapes of building blocks, the electrostatic interactions between them, and the morphology of the resulting structures. 25 For example, proteins having different amino acid sequences can fold into very similar shapes, and subsequently self-assemble into oligomers and other hierarchical structures, such as fibres, closed shells, or tubes. 26,27 . Further examples are the multicellular tumor spheroid (MCTS) models for mimicking the microenvironment of tissues. 28 These experiments have shown the effect of surface charge on nanoparticle penetration into a MCTS.
Directed self-assembly of polarizable ellipsoids in an external electric field has been computationally studied using Monte Carlo simulations of a two point-charge model of polarizable prolate ellipsoids. 29 However, there have not been corresponding developments towards a general methodology for treating electrostatic interactions between non-spherical particles. Exact solutions to this problem have only been presented for a single uniformly charged spheroidal shell, 30 and where the image charge method has been used to treat conducting ellipsoidal particles. 31 In this paper, an analytical theory of electrostatic interactions between spheroidal particles has been developed, building on previous work, [32][33][34][35][36] where analytical expressions have been given for the electrostatic force between charged, dielectric sphere -sphere 32 and sphereplanar surface systems. 33 In these electrostatic models, the mutual effect of charge is obtained from Gauss's law, which couples uniquely the electrostatic potential with the distribution and magnitude of electric charge on the surfaces of the interacting objects. The accumulated sur-face charge is integrated to obtain an analytical expression for the electrostatic force acting on interacting objects at arbitrary separation. The result is a simple series expression for the force that can be efficiently generalized for studying interactions not only in vacuum 32,33 but also in solution 34 and in electrolytes. 35,36 The solution has been evaluated by comparison with existing solutions for a range of simple geometries including a point charge corresponding to a non-polarizable sphere, a charged rod corresponding to a non-polarizable prolate spheroid and a disc corresponding to a non-polarizable oblate spheroid.

II. METHODOLOGY
A. Geometry the problem and expansion of the electrostatic potential The problem to be addressed involves two dielectric spheroidal particles, denoted as i = 1, 2 in Fig. 1, of arbitrary size and defined by semi-axes a i and c i , permittivity k i , and carrying an arbitrary charge Q i in a surrounding dielectric medium of permittivity k m . The particles are placed on the same axis of symmetry z at the distance R between their centres. The problem is solved in spherical coordinate systems with an origin at the centre of the spheroids. The distribution of electric potential inside and outside the spheroids is described by the Laplace equation which is supplemented by two boundary conditions. The first assumes continuity of the electric potential on the surface of the i-th spheroid: where Φ i,in is the potential inside the spheroid, and Φ i,out + Φ j,out is the potential outside the spheroid with contributions from both the i-th and j-th spheroids, j = 3 − i, r i is the radial coordinate in the spherical coordinate system with the pole in the center of the i-th particle, ρ i (µ i ) is the spheroid surface radial coordinate in the spherical frame system: The second boundary condition states that the normal component of the dielectric displacement field is discontinuous due to the presence of a free charge on the surface of a spheroid: Here n i is the unit normal vector on the surface of the i-th spheroid, σ i (µ i ) is the surface charge density of the i-th spheroid, and ε 0 is the permittivity of vacuum.
The electrostatic potential inside the i-th spheroid, which satisfies the Laplace equation (1) can be expanded in terms of Legendre polynomials P n (µ i ): 37 The potential outside each spheroid that satisfies equation (1) and vanishes at infinity, takes the form 37 In order to apply boundary conditions (2) and (3) and determine the expansion coefficients A n,i and B n,i it is necessary to re-expand the potential (5) and use only one set of spherical coordinates for each spheroid 38 The corresponding derivatives of the electrostatic potential are B. The case of isolated spheroid The surface charge distribution σ i (µ i ) is found from the assumption that the surface of an isolated spheroid is equipotential: where the surface potential φ 0 is described as: 39,40 FIG. 1: A geometric representation of two interacting dissimilar spheroids. Dielectric constants, permanent charges, and the semi-axes for spheroids 1 and 2 are denoted as k 1 , Q 1 , a 1 , c 1 and k 2 , Q 2 , a 2 , c 2 .
Using expansion (5) in (8) gives: where the B n are constant coefficients corresponding to an isolated spheroid. Expanding both parts of equation (10) in terms of Legendre polynomials yields where Solution of the linear system (11) gives the expansion coefficients B n,i . Inside an isolated spheroid with an equipotential surface the electric field is zero; therefore the second boundary condition, (3), can be rewritten as Using (5) and the expansion coefficients B n,i , the surface charge distribution is given by: Here and thereafter, the components of the normal vector n = n rr + n θ 1 − µ 2θ on the surface of the spheroid are defined as Note that equation (14) is expressed in a general form and can be applied to any three dimensional shape with axial symmetry. In this paper, it is tested against the known formula for the surface charge density on an isolated spheroid with a uniformly distributed potential 39 Figure 2 compares the numerical results obtained using Eq. 14) and the analytical expression (16) for three different cases corresponding to the aspect ratio of a : c = 1 : 1 (sphere), a : c = 3 : 4 (prolate), and a : c = 4 : 3 (oblate). The deviation of the numerical results is within 0.1%, mainly in the charge deficient areas as compared to the distribution of charge on the surface of a sphere; thus demonstrating the reliability of the proposed method.

C. Two spheroids at a finite separation
If two spheroids are located at a finite distance apart, the boundary condition (2) takes the following form: Here, the electrostatic potential of the j-th spheroid is re-expanded in a spherical coordinate system with the origin at its centre using an addition theorem for Legendre polynomials. 38 Multiplying both sides of (17) by P k (µ i ) and integrating over the limits −1 to 1 yields: where and β kn,i is defined by Eq. (12). The second boundary condition (3) expanded in terms of Legendre polynomials takes the form: where  14) are shown by circles (sphere), triangles (prolate) and squares (oblate). The embedded plots represent the relative deviation between the two approaches. and Substitution of the derivatives (7) of the electrostatic potential into Eq. (21) gives: where the following notation has been introduced: Similarly, substitution of Eq. (14) into Eq. (22) yields: Hence, Eq. (20) can be rewritten as Finally, combining Eqs. (18) and (26) gives the required set of linear equations for the coefficients A n,i and B n,i : where i = 1, 2, j = 3 − i, k = 0, 1, 2, . . . , ∞,

D. Electrostatic force
The Maxwell stress tensor is used to calculate the electrostatic interaction force 41 acting on spheroid i due to the presence of spheroid j where is the normal component of the Maxwell stress tensor and τ = −n θ 1 − µ 2 e r + n r e θ is the tangent unit vector.
Here, the normal component of the electric field is given by: The tangential component of the electric field is defined as: Eq. (30) can be substituted into Eq. (29) for the electrostatic force and rewritten as: which solves the posed electrostatic problem.

III. RESULTS AND DISCUSSION
We next consider the effect of non-sphericity on the nature of electrostatic interactions between two polarisable spheroids of the same shape, size (a 1 = a 2 ≡ a, c 1 = c 2 ≡ c) and dielectric constant (k 1 = k 2 ≡ k), but with different charges Q 1 /Q 2 = 2, whilst keeping the capacitance of spheroids constant. We assume that at infinite separation distance between spheroids their capacitance is equal to the capacitance of a sphere. The assumption of the constant capacitance implies that during the deformation of an isolated sphere the ratio between its surface charge and surface potential remains constant. Therefore, this assumption has been chosen as the most physically meaningful for the case when the effect of non-sphericity on the electrostatic interaction is studied and the effects of changed charge and/or potential are excluded.
This approach allows us to find the relationship between the axes of a spheroid, a and c, and radius of the corresponding sphere, r, into which the spheroid degenerates at a = c. Under the assumption of constant capacitance where and Using the non-sphericity parameter x = a/c and substituting Eqs. (35) and (36) into (34) gives the following relationship between the radii of a spheroid, a, c and the radius of the equivalent sphere, r: The relationship (37) has been used in all numerical tests presented in this paper for a range of x from 0.83 to 1.17. The electric potential generated by a point charge Q point is typically represented by equipotential surfaces (regions in which every point has the same potential), which take the form of concentric spheres centred at the point charge. 39 If the point charge is substituted by a small sphere with the same charge, Q sphere = Q point , and the uniformly distributed surface potential (this condition also implies the uniformly distributed surface charge), the sphere will create the same electric potential outside its boundaries as the point charge (Gauss's law). Therefore, the electrostatic forces between two point charges and two uniformly charged non-polarizable and non-overlapping spheres are equivalent. The same reasoning can also be applied to charged objects of any arbitrary shape using the superposition principle. For example, a prolate spheroid has an equipotential surface of a uniformly charged rod with the length equal to the interfocal distance 2f of the corresponding spheroid, and an oblate spheroid has an equipotential surface of a disc with the radius equal to the radius of the focal line f and with the following radial distribution of surface charge density: 41 where Q disc is the charge of the disc, and r d is a radial coordinate on the disc surface. Table I contains the analytical equations of the electrostatic forces for the following cases: • two point charges at a distance R, which corresponds to two spheres of radii r 1 and r 2 at a surfaceto surface separation s = R − r 1 − r 2 ; • a charged rod of length 2f and a point charge separated at a distance R from the center of the rod that corresponds to a prolate spheroid with c = f / √ 1 − x 2 and a sphere of radius r at a surfaceto surface separation s = R − c − r; • a charged disc of radius f and a point charge at a distance R from the center of the disc that corresponds to an oblate spheroid with c = f / √ x 2 − 1 and a sphere of radius r at a surface-to surface separation s = R − c − r; • a charged disc of radius f 1 and a charged rod of length 2f 2 at a distance R between their centers that corresponds to an oblate spheroid with c 1 = f 1 / x 2 1 − 1 and a prolate spheroid with 2 at a surface-to surface separation s = R − c 1 − c 2 ; • two charged rods of lengths 2f 1 and 2f 2 at a distance R between their centers that corresponds to two prolate spheroids with c 1,2 = f 1,2 / 1 − x 2 1,2 at a surface-to surface separation s = R − c 1 − c 2 .
The electrostatic force calculated using Eq. (33) for nonpolarizable spheroids: k 1 = k 2 = k m , including a sphere as a specific case, should, therefore, give the same result as the electrostatic force obtained from the simple expressions summarised in Table I for cases involving a charged rod, a disc, and a point charge. For these simple geometries, Fig. 3 compares calculations of the electrostatic force as a function of surface-to-surface as defined with reference to Table I. The results obtained using the methodology presented above and the corresponding analytical expressions given in Table I are in excellent agreement. These limiting cases can be interpreted as electrostatic forces between non-polarizable spheroids. Consider next the transition from repulsion to attraction of like charged identical spheroids by changing their eccentricity. It should be noted that determining the range of parameters in which, for particles of the same charge, a transition from repulsion to attraction takes place as the distance between them decreases is not a trivial task. Even for spherical particles, the boundaries of this region are determined by the ratio of the charges on the particles, the ratio of their sizes and their dielectric constant relative to that of the medium. 32,42 For particles of a spheroidal shape, even for the case of an axially symmetric distribution of surface charges, the eccentricities have to be added to the parameter space. This challenging task remains outside the scope of the present work, which is focused primarily on the development of analytical and numerical solutions to the electrostatic interactions between different non-spherical geometries of charged dielectric particles.
It is well known that for the case of spherical particles of the same radius carrying equal charges, there is no attraction even for conducting particles. 42,43 Consider the case when one particle carries twice the charge of the other, i.e. Q 1 = 2Q 2 . As the distance between the surfaces of the particles decreases, the number of terms in the multipolar expansion required for an accurate estimation of the potential increases and the dimensionality of the set of algebraic equations defining the expansion coefficients increases accordingly (see 44 ). Therefore, the test calculations are restricted to interactions at sufficiently large interparticle distances, s = 0.01r, at which for spherical particles with the same dielectric constant repulsion transforms into attraction at k 1 = k 2 ≈ 18.5 (k m = 1). Therefore, for comparison we consider the values k 1 = k 2 = 18, 18.5 and 19.  . 3: (a) Electrostatic force, scaled by the Coulomb force, calculated as a function of the surface-to-surface separation for a prolate spheroid and a sphere (1 ), an oblate spheroid and a sphere (2 ), an oblate and a prolate spheroids (3 ), two prolate spheroids (4 ), and two spheres (Coulomb force) (5 ). Lines are analytical results given by the equations in Table I, symbols are numerical calculations using Eq. (33) with x = 4/3 for oblate spheroids and x = 3/4 for prolate spheroids; (b) Relative errors.
For the case of polarizable spheroidal particles in vacuum (k m = 1), Fig. 4 shows the electrostatic force between two identical spheroids carrying different amount of charge (Q 1 /Q 2 = 2) as a function of the non-sphericity parameter, x, calculated for three values of the dielectric constant, k = 18, 18.5 and 19. The spheroids are kept at a fixed surface-to-surface separation s = 0.01r. The values given for the dielectric constants have been selected from extensive numerical experiments, to reveal the switch in electrostatic behaviour from attraction to repulsion, between like-charged spheroids depending on the value of the non-sphericity parameter, x. For k = 18.5 and above, the interaction can switch from a counter-intuitive attraction between like-charged particles (negative value of the force ratio) driven by charge-induced polarisation, to repulsion (positive value of the force ratio). This switch occurs either as the shape of the interacting spheroidal particles changes from oblate to prolate or if the value of the non-sphericity parameter for two oblate spheroids is increased sufficiently. For chosen values of dielectric constant and charge ratio, the minimum in the electrostatic force corresponds to two oblate spheroids (x > 1). This behavior is as a result of a specific distribution of the surface charge, which depends on the non-sphericity parameter (see Fig. 2), and the complex nature of polarisation interactions between spheroids.

IV. CONCLUSIONS
An analytical expression for the electrostatic force acting between two dielectric spheroids located on the same axis of symmetry is presented. Variation in the electrostatic force with a change in the value of the nonsphericity parameter shows an interesting switch in electrostatic behaviour between two like-charged spheroids with a charge ratio of 2. At a critical value of the dielectric constant, k = 18.5, and above, the F/F Coulomb ratio has a negative value, which corresponds to an attractive interaction between like-charged oblate spheroids. If the shape of the interacting spheroidal particles changes from oblate to prolate or if the value of the non-sphericity parameter of two oblates is increased sufficiently the interaction switches from attraction to repulsion. The proposed analytical model was benchmarked against existing analytical solutions for the interaction between non-polarisable rods, discs, and point charges, and against an earlier electrostatic model for dielectric shperes 32 showing excellent agreement. The result is of practical significance and represents a first step towards a more general theory of electrostatic interactions between non-spherical objects as it can be generalized to any arbitrary shape with axial symmetry, as shown in the approach taken in Ref. 45 . Derivations for the electrostatic force for the simple limiting cases of a charged rod, disc, and point charge and additional computational issues are discussed in Appendices.

V. CONFLICTS OF INTEREST
There are no conflicts to declare.

ACKNOWLEDGMENTS
Finacial support from a ERC Consolidator grant (ERC-2012-StG 307755 FIN), the Innovate UK grant (102696) and the International Space Science Institute (Bern) are gratefully acknowledged. A. J. S would like to thank the Leverhulme Trust for the award of an Emeritus Fellowship.
Appendix A: Derivation of electrostatic forces for simple limiting cases given in Table I 1. Uniformly charged rod and point charge The force between a uniformly charged rod and a point charge (see Fig. 5) can be derived by integrating the force between an infinitely small element of the rod and the point charge over the length of the rod: where K = 1 4π 0 ≈ 9 × 10 9 Vm/C is a constant of proportionality.
FIG. 5: A uniformly charged rod with charge Q rod and length 2f and a point charge Q point . R is the distance between the centre of the rod and the point charge.

Two uniformly charged rods
The force between two uniformly charged rods (see Fig. 6) can be derived by a double integration of the force between infinitely small elements of the rods over their lengths: Here z = 0 is assumed to be the centre of the first rod.

Charged disc and point charge
The force between a uniformly charged ring and a point charge (see Fig. 7) can be derived by integrating the force FIG. 6: Two uniformly charged rods with charges Q rod1 and Q rod2 and lengths 2f 1 and 2f 2 . R is the distance between their centres.
between an infinitely small element of the ring and the point charge over the circumference of the ring: FIG. 7: A uniformly charged ring with charge Q ring and radius r and a point charge Q point . R is the distance between the centre of the ring and the point charge.
The force between a disc with the surface charge density and a point charge (see Fig. 8) can be calculated by integrating the force from Eq. (A3) between an infinitely thin ring element of the disc and the point charge over the circumference of the ring:

Charged disc and uniformly charged rod
The force between a disc with the surface charge density and a uniformly charged rod (see Fig. 9) can be calculated by integrating the force from Eq. (A4) between the FIG. 8: A charged disc Q disc with the radius f and a point charge Q point . R is the distance between the centre of the disc and the point charge.
charged disc and an infinitely small element of the rod over the length of the rod: FIG. 9: A charged disc Q disc with the radius f 1 and a uniformly charged rod Q rod with the length 2f 2 . R is the distance between the centres of the disc and the rod.
Appendix B: Limiting case of two polarizable spheres As verification, the presented methodology has been tested for the case of two dielectric spheres: a i = c i ≡ a i , i = 1, 2. In this case φ i,0 = Q i 4πk m ε 0 a i and ρ i (µ i ) = a i . Eq. (11) takes the form: which gives The normal vector components (15) are deduced to Substituting (B2) and (B3) to (14) gives: Therefore, the first equation in (27) takes the form: whereas the second equation in (27) gives: Substituting (B5) in to (B6) yields: The same equation was obtained in Ref. 32 for the case of k m = 1.
In order to confirm the analytical derivations made above, the electrostatic force between two polarizable spheres of the same size has been calculated using Eq. (33) and compared with the force calculated using the model from Ref. 32 . The following parameters have been chosen to reproduce the most relevant case of like charge attraction: k 1 = k 2 = 20, Q 1 = 1e, Q 2 = 10e, and k m = 1. Figure 10 demonstrates that the results are in a good agreement and the relative error does not exceed 1% at short separation and less than 10 −5 % when the spheres are far apart. The greatest error (0.15%) is achieved at the point where the force changes sign and crosses the x-axis.
Appendix C: Error analysis related to the number of terms in the expansion of the electrostatic force Convergence of the present methodology is demonstrated for the example of two geometrically identical spheroids with the same dielectric constants but different charges. Numerical experiments showed that in the case of Q 1 /Q 2 = 2 and k 1 = k 2 = 17 and 20 the solution is stable for 0.83 ≤ x ≤ 1.17. Figure 11 shows the electrostatic forces and the calculation errors versus number of terms for spheroids described here at a separation of 0.01r. For the examples of prolate spheroids and spheres the method shows excellent convergence at values of n in the range 20 ≤ n ≤ 140, whereas for the case of oblate spheroids convergence stops at n = 120 and thereafter the error increases. The linear system (27) is generally sparse and ill-conditioned, i.e. contains many zero elements and elements with large differences in values. Moreover, the problem of two oblate spheroids has no trivial solution FIG. 10: Electrostatic force (relative to the Coulomb force) between two polarizable spheres of the same radius r with the dielectric constants k 1 = k 2 = 20 and charges Q 2 /Q 1 = 2 in vacuum k m = 1 calculated by means of the methodology from Ref. 32 (line) and the present model (33) (symbols) versus surface-to-surface separation s relative to the sphere radius r. The embedded plot represents the relative difference between the forces.
for the non-polarizable case (unlike the problem of two prolate, non-polarizable spheroids). Three separate numerical methods have been tried to solve the problem (27): 46 lower-upper (LU) decomposition with iterative improvement of a solution, singular value decomposition (SVD) and preconditioned biconjugate gradient method (PBCG). All methods give identical results and for further calculations the LU-decomposition method has been chosen.