Screened Coulomb interactions of general macroions with nonzero particle volume

A semianalytical approach is developed to calculate the effective pair potential of rigid arbitrarily shaped macroions with a nonvanishing particle volume, valid within linear screening theory and the mean-ﬁeld approximation. The essential ingredient for this framework is a mapping of the particle to a singular charge distribution with adjustable effective charge and shape parameters determined by the particle surface electrostatic potential. For charged spheres, this method reproduces the well-known Derjaguin-Landau-Verwey-Overbeek (DLVO) potential. Further exemplary benchmarks of the method for more complicated cases, like tori, triaxial ellipsoids, and additive torus-sphere mixtures, leads to accurate closed-form integral expressions for all particle separations and orientations. The ﬁndings are relevant for determining the phase behavior of macroions with experiments and simulations for various particle shapes.


I. INTRODUCTION
Screened Coulomb interactions of electronic or ionic nature are ubiquitous in quantum-mechanical and classical systems, such as strongly correlated electron matter [1,2], chemical bonds [3,4], superconductors [5], proteins [6], liquid crystals [7][8][9], DNA [10][11][12], graphene [13], lipid membranes [14,15], supercapacitors [16,17], microfluidics [18], and dusty plasmas [19,20].A general understanding of electrostatic screening in various geometrical settings is needed, considering that many of these systems have a complex geometry.In particular, classical charge-screened particles of various shapes and surface functionalities can be experimentally synthesized and characterized today in great detail [21,22]; however, theoretical understanding of effective particle interactions is lagging behind as it is difficult to account for finite particle volume and nonspherical particle shape.This paper is aimed at bridging the gap between the available experimental and theoretical toolkits.
In order to express the system in solely the degrees of freedom of interest-such as the positions and orientations of specific particles-it is useful to integrate out the "fast" charge degrees of freedom, which leads to an effective description in terms of electrostatic screening.In free-electron like metals, this procedure leads to Friedel oscillations, electron density modulations near a solid-fluid interface and around impurities [23][24][25], which is a canonical example of an emergent phenomenon caused by screening.In classical systems [26] and sufficiently dilute quantum systems [27,28], screening is often associated with the damped spatial decay of the electrostatic potential (and thus the effective pair potential), that for point particles and for low voltages compared to the thermal energy has the Yukawa form, ∼ exp(−κr)/r, with κ −1 being the Debye screening length and r being the radial distance, as opposed to the bare Coulomb case ∼1/r.Of special interest are spherical charged particles, with radius a in the colloidal (sub)micron regime, dispersed in ion-containing liquids, because of their tunable charge and screening properties [29].Integrating out the degrees of freedom of the smaller ions results in a Yukawa-type effective sphere-sphere potential with a prefactor that depends on the particle charge and, unlike for point particles, also the salt concentration via κa arising from the ion-impenetrable particle hard core [30,31].This so-called Derjaguin-Landau-Verwey-Overbeek (DLVO) potential [32] is an essential theoretical tool for understanding the behavior of charge-stabilized colloids [33], even for out-of-equilibrium suspensions [34].
For nonspherical shapes, the screened-electrostatic pair interaction is only analytically known in a few cases for all particle configurations, even within linear screening theory.However, some studies exist for disks [35][36][37][38], rods [39,40], spheroids [41][42][43][44][45], or helices [46], where the potential is sometimes calculated only for infinitely long, thin, or ion-penetrable particles, restricted particle configurations, or orientation-averaged interactions [47].The difficulty in finding analytical solutions lies in the finite ion-impenetrable particle volume which complicates matching the series expansion solution (if it is even available for the geometry under consideration) of the unscreened potential inside and the screened potential outside the particle via the boundary conditions.However, when the pair potential would be known, one does not need to numerically solve the three-dimensional Poisson(-Boltzmann) equation for every single-particle configuration at each simulation step, as in Refs.[48][49][50] for simulations of charged colloids.Instead, computationally less expensive simulations with effective pair potentials can be used, and by mapping to cell models, even charge regulation can be incorporated [51].The lack of availability of accurate pair FIG. 1. Scheme of mapping an arbitrary ion-impenetrable particle with surface charge density σ to an effective ion-penetrable charge distribution q(r), that can be a point, line, or surface charge.potentials might explain why fewer phase behavior studies are known for complex-shaped charged particles [52][53][54] than for charge-neutral hard particles [55][56][57][58][59][60].
In this paper, I devise a framework to semianalytically approximate effective interactions between (not necessarily equal) charged finite-size particles with not necessarily spherical shape.By mapping particles to singular charge distributions (i.e., expressed by δ functions), I find a straightforward and accurate evaluation of the interaction free energy for arbitrary interparticle separations and orientations.For spheres, this method reproduces DLVO theory, and hence a similar level of approximation is expected such as weak double layer overlap [61].After discussing spheres, I apply the framework to more complicated shapes, such as tori and ellipsoids.

II. GENERAL THEORETICAL FRAMEWORK
To set up the theoretical framework, I consider an ionimpenetrable charged particle of arbitrary shape with dielectric constant p , surface P, and interior volume int(P ), immersed in a structureless solvent with dielectric constant s and Bjerrum length B = βe 2 /(4π 0 s ), with e being the proton charge, 0 being the vacuum permittivity, and β −1 = k B T , with k B being the Boltzmann constant and T being temperature; see the scheme in Fig. 1.For simplicity, I focus on a 1 : 1 salt, such that the mean-field approximation is valid.By thermally averaging over the ions, an inhomogeneous electrostatic potential φ(r)/(βe) describes the electric double layer that is formed around the particle.I split the total electrostatic potential as a contribution inside the particle φ < (r) = φ(r)| r∈int(P ) , and a contribution outside the particle φ > (r) = φ(r)| r / ∈int(P ) .For linear screening, there is the condition on the dimensionless electrostatic potential |φ(r)| 1, meaning that I only consider potentials much smaller than the thermal voltage.Within the mean-field approximation, φ < (r) and φ > (r) are given by the Laplace and Debye-Hückel (DH) equations, respectively, with κ −1 = (8π B ρ s ) −1/2 being the Debye screening length where ρ s is the reservoir salt concentration, and I enforce continuity φ < (r) = φ > (r) for r ∈ P, and constant-charge boundary condition on the particle surface with outward normal ν and surface charge density eσ , The finite particle volume complicates matching φ > (r) and φ < (r) through the boundary conditions [42], although formal exact multiple-scattering expansions exist [62].Sometimes Eq. ( 1) is not even separable in certain coordinate systems, which further complicates finding analytical solutions.For example, due to the nonseparability of the Helmholtz equation (DH with imaginary κ) in toroidal coordinates, the complete solution can only be expressed in terms of lengthy toroidal wave functions [63].However, some approximations for the double layer around a torus exist [64,65].
The first central message of this paper is that an ionimpenetrable charged particle can be mapped to a singular charge distribution q(r) without a particle hard core, described by the DH equation (∇ 2 − κ 2 )ϕ(r) = −4π B q(r), such that the outside potential is approximated to a very high accuracy by φ > (r) ≈ ϕ(r), schematically shown in Fig. 1.Consequently, ϕ(r) can be expressed as the convolution of q(r) with the DH Green's function, Here, q(r) can either be a point, line, or surface charge distribution and has to be parametrized with the same symmetry of the particle.The salt-dependent shape and magnitude of q(r) can be determined by matching, e.g., the numerically or analytically obtained surface potential.The benefit of this method is that the same q(r) for a single particle enters the approximate analytical expression for the charge-screened two-particle interaction, which is the second central message of this work.For two arbitrary particles with center-to-center distance vector d and orientations 1,2 that are mapped to charge distributions q i (r; d, i ) (i = 1, 2), respectively, I find for the effective pair interaction within the DH approximation where I subtracted the infinite self-energy terms, and for simplicity I omitted the explicit configurational dependence of the charge distributions in the right-hand side of the equation.

A. Charged spheres
As an example, a sphere of radius a and total charge Z s e = 4π a 2 σ e can be mapped to a point charge q(r) = Q p δ(r).Using Eq. (3) gives ϕ(r) = Q p B exp(−κr)/r, to be compared with the exact analytical solution [66] Matching the potential on the particle surface, ϕ(r = a) = φ > (r = a), I find Q p = Z s ϒ s , with ϒ s = exp(κa)/(1 + κa).Alternatively, Q p can be computed: For the ion densities caused by the singular point charge ρ ± (r), one can check that showing that an ionimpenetrable charged sphere produces the same electrostatic potential for r > a as an ion-penetrable particle, consisting of a suitable point charge surrounded by a plasma of ions.This physical interpretation is possible because φ < (r) = Z s B /[a(1 + κa)] is constant and therefore does not contribute to Eq. ( 2).
For spherical particles, the DLVO potential results from Eq. ( 4), with d = |d|.The same result follows from the linear superposition approximation (LSA) used on Eq. ( 5), φ 2B (r) ≈ i φ > (r − X i ), in force calculations with the stress tensor [66].However, the free energy route as in Eq. ( 4) is incompatible with the LSA: One would obtain a different result than Eq. ( 6) because the LSA does not account for ion-particle hard-core interactions [36].However, unique for the pointcharge mapping is the compatibility with the free energy and the stress tensor route, because the particle hard core is effectively mapped out.Finally, the choice of q(r) is not unique: Mapping spheres to ion-penetrable charged spherical shells gives the same result as the point-charge mapping.

B. Charged tori
Next, I show an example where the particle is mapped to a line charge C. Such a mapping is akin to the slender-body theory for creeping flow [67]; however, my method works also for "thick" particles, because C need not coincide with the centerline of the particle.As an illustration, I consider the pair potential between two identical tori with inner radius R I , outer radius R O [Fig.2(a)], and uniform surface charge density σ = Z t /(4π 2 R I R O ).As a first step, I map the torus to a charged ring C with parametrization γ : [0, 2π ) → C given by γ (u) = (R cos u, R sin u, 0), of uniform line charge density λ = Q r /(2π R).Unlike the point-charge mapping, not only does the line charge number Q r have to be determined, but also the shape for R ∈ where I factorized the ring charge Q r = Z t ϒ t using the linearity of Eq. ( 1).Note that only for κ = 0 (and along the z axis) can the integral in Eq. ( 7) be evaluated in terms of known special functions [68].For the general case, the integral can be numerically computed on a desktop PC within seconds.Second, to establish the values of ϒ t and R, I fit the surface po- intraparticle double-layer overlap, Eq. ( 7) agrees excellently with numerics, capturing the full spatial dependence of φ > (r) for suitably chosen ϒ and R, even for tori with large particle volumes.Finally, Eq. ( 7) gives analytical insight: For r → ∞, I find that ϕ(r) ∼ A(θ )Z t B exp(−κr)/r, with anisotropy function (see Appendix D for the derivation), showing the well-known result that particle anisotropy still persists in the far-field electrostatic potential unlike the unscreened case, and that the anisotropies are more pronounced for large κ [36,42].To my best knowledge, A(θ ) has never been calculated for a torus before.
To calculate effective pair interactions, I parametrize two identical rings with arbitrary orientations as and ni defined in the inset of Fig. 3(b).Using Eq. ( 4), I find the interaction between two identical tori Note that the difference |γ 1 (u) − γ 2 (v)| depends on d, n1 , and n2 .Despite being an integral representation, the evaluation of Eq. ( 9) is far less computationally expensive than threedimensional finite-element calculations of Eq. (1) for every fixed particle configuration.Explicit calculations give for d not parallel to n1 and/or n2 , where the term f in the integrand simplifies to whereas for d not parallel to either n1 or n2 , In Eq. ( 12), I used the parametrization ˆli In Fig. 3, I compare the numerical calculation and the analytical approximation of Eqs.(10) for a wide variety of particle configurations and Debye screening lengths.I find an excellent agreement with only a deviation for particle separations close to contact where the repulsion is underestimated, as expected from a method that is equivalent to the LSA [61].

C. Additive torus-sphere mixtures
The method can also be extended to additive mixtures of particles, as I highlight for the torus-sphere interaction which, again, agrees excellently with numerics (Fig. 4), even when the sphere partly enters the hole of the torus [Figs.4(b) and 4(c), blue line].In the insets, I highlight on a log-linear scale that the decay length for these particle shapes is still κ −1 ; however, with an orientation-dependent interaction amplitude of higher anisotropy when the salt concentration is increased [Figs.4(b) and 4(c)], which is generic for anisotropic particles [36,42], as is for the electrostatic potential [Eq. ( 8)].

D. Charged triaxial ellipsoids
To show the generality of the framework, I consider the interaction potentials of two triaxial ellipsoids, for which a mapping to a spherical shell is needed, rather than a line charge.I hypothesize, however, that some metal spheroidal particles can be mapped to straight lines, because the (unscreened) isopotential surfaces of straight lines are prolate spheroids with specific aspect ratios [70].It is straightforward to generalize the approach to mappings to ion-penetrable charged surfaces S with surface charge distribution σ s , and S parametrized by : [a, b] × [c, d] → S. I find from Eq. ( 5) where g(u, v) is the determinant of the induced metric tensor.
For the particle, consider the parametrization for a triaxial ellipsoid X(φ, θ ) = (a l cos φ sin θ, a m sin φ sin θ, a n cos θ ), (15) where the angles are defined for φ ∈ [0, 2π ) and θ ∈ [0, π ), with volume V e = (4/3)π a l a m a n , and surface area with F (ψ, k) and E (ψ, k) being the incomplete elliptic integrals of the first and second kinds, respectively.Furthermore, cos ψ = a n /a l and Since X and should have the same shape, I use an ellipsoidal shell with the same aspect ratios as the particle, but with a smaller size parameterized by R l , with u ∈ [0, 2π ), v ∈ [0, π ), γ m = a m /a l , and γ n = a n /a l .From Eq. ( 14), I determine the electrostatic potential of a triaxial ellipsoidal particle where I used for σ s a homogeneously charged ellipsoidal shell with total charge Z e ϒe.The parameters R l and ϒ are determined via a fit of the surface potential ϕ 0 (φ, θ ) = ϕ(X(φ, θ )).As an example, I determine ϒ and R l for a triaxial ellipsoid with the  18)] and (d) the effective pair potential [Eq. ( 21)] of the semianalytical ellisoidal-shell mapping with numerical finite-element calculations of the Poisson-Boltzmann equation for a triaxial ellipsoidal particle at various values of κa l , with a l being the distance from the origin to the particle surface along the x axis.The potential is compared along various cuts around the ellipsoid, as well as an enlargement of the comparison between numerical and analytical results close to the particle surface.In all plots, p / s = 0.2, B /a l = 0.01, and Z = 100.The shell that has been used to represent the particle has the same ratios of the ellipsoidal main axes, but with a rescalable parameter R l < a l and charge parameter ϒ with the same meaning as for the torus.For precise dimensions of the particle, see the main text.
same volume and surface area as a torus with inner radius R I = 0.5R O and a l = R O .In this case, the aspect ratios are γ m = 2.595 and γ n = 0.454.Moreover, I fix κa l = 1 and p / s = 0.2.The result of the fit is R l = 0.9599 and ϒ = 1.1523.In Figs.5(a) and 5(c), I show a comparison of the resulting electrostatic potential with finite-element calculations, showing excellent agreement outside the particle.
Finally, by using Eq. ( 4) I derive for the effective pair interaction between two particles that are mapped to charged shells, I apply Eq. ( 19) to arbitrarily oriented and positioned ellipsoidal particles, mapped to (for i = 1, 2) Here { ˆli , mi , ni } is an orthonormal triad of vectors.I find applying by Eq. ( 20) to Eq. ( 19) that where the function f in the integrand is expressed as with ˆli = mi × ni .This ellipsoid-ellipsoid interaction potential is tested in Fig. 5(d) to numerical calculations, showing excellent agreement, except for small particle separations, as was also found in the previous examples.

IV. CONCLUSIONS AND OUTLOOK
In summary, I developed a framework to derive effective pair potentials between finite-sized arbitrarily shaped rigid macroions.With different combinations of spheres, ellipsoids, and tori in various mutual orientations, I showed the applicability and accuracy of this method: This framework gives analytical insights with broad relevance for experiments and simulations.I expect that the method applies to many shapes when multiple inhomogeneous singular charge distributions are used.Furthermore, as recently shown, the method works also for spheres dispersed in nematic liquid crystals [9].Finding the correct charge distribution for a specific particle can be nontrivial, but often symmetry arguments and an analysis of the necessary multipole moments to match the far-field electrostatic potential are useful considerations.For example, it was suggested that Janus spheres can be described by a collection of point charges [71].My method becomes less accurate at small particle separations, like DLVO theory, where the surface potential and/or surface charge density become "polarized," which, however, can be reconciled with the Derjaguin approximation [43,66,72,73].Furthermore, it is not possible to derive multibody interactions (of the form as in Ref. [74]), because of the underlying LSA equivalence.
As an outlook, I propose extending the theory with (manybody) charge regulation [51,73,[75][76][77] and renormalization [78][79][80][81] to incorporate more types of electrostatic boundary conditions and nonlinear screening, respectively.The expressions of this paper can then still be used with the bare charge replaced by an effective (renormalized) charge.Another extension of the theory would be to include correlations beyond the mean-field result presented here.It would be interesting to see in this case whether the mapping to a singular charge distribution is still suitable to describe the (thermally averaged) electrostatic potential and pair interactions accurately.Such questions can be answered via the field-theoretical formulation [82] of two particles immersed in an ion-containing solvent.Furthermore, such a fieldtheoretical formulation would be interesting to investigate the effects of flexible particles, rather than the rigid particles discussed here.In this case, an elastic free-energy contribution would self-consistently determine the singular charge distribution needed for the mapping, and this deserves further research.
More broadly, the findings might also be valuable for any physical system governed by the Helmholtz equation, e.g., acoustics [83] and optics [84], or systems where the Yukawa potential is involved, such as wetting [85].Finally, it would be intriguing to explore charge-screened active matter for various "thick" particles.Here, the Yukawa potential is already often used to model steric repulsions between thin active rods [86,87].FIG. 6. Examples of fits of the numerically obtained surface potential (scaled with the Bjerrum length B and particle charge Z) with the expression of a ring charge mapping, Eq. (B2).From this fit the values for R and ϒ are obtained.All calculations were performed at To numerically evaluate the expressions obtained from a singular-charge distribution mapping, such as Eqs.( 9), (11), and ( 14) from the main text, and Eq. ( 18) in this document, I used the standard integral and integral2 commands in MATLAB.For the ellipsoid-ellipsoid interaction Eq. ( 21), iterative use of these commands was made.

APPENDIX B: DETERMINATION OF CHARGE AND SHAPE PARAMETERS FOR AN ELECTROSTATICALLY SCREENED CHARGED TORUS
Recall the approximation of the electrostatic potential for a single torus Eq. ( 7) written down in cylindrical coordinates, which for suitably chosen values of R and ϒ approximately describes a toroidal particle P with standard parametrization ).Consequently, I find the approximation for the surface potential In MATLAB, I use the fit command from the Curve Fitting Toolbox to establish the values of ϒ and R by fitting Eq. (B2) to numerically obtained surface potentials from COMSOL using the Levenberg-Marquadt algorithm.A few examples of fits are shown in Fig. 6.The fits are of good quality even for large particle volumes.The quality of the approximation, however, seems to be the least accurate for extremely thin double layers κR O = 10, although this turns out to be of little importance for the accuracy of the full spatial dependence of the electrostatic potential around the particle, as I shall demonstrate later.It is meaningful to explore the values of the charge parameter ϒ and shape parameter R for a wide variety of particle thicknesses as quantified by R I /R O and screening lengths as quantified by κR O .For this I fitted a whole set of numerically obtained surface potentials to Eq. (B2).Regarding ϒ, it is meaningful for the comparison between different κR O and R I /R O to scale out a large factor for better comparison, Result of (a) the charge parameter and (b) shape parameter R as obtained from fits of the surface potential, see Fig. 6, for various values of κR O and torus thicknesses as quantified by R I /R O .The results are independent of particle charge Z and Bjerrum length B in the linear screening regime, but depend on p / s (fixed at 0.2 as in Fig. 6).In panel (c), we plot a statistical parameter showing the quality of the fits, the so-called coefficient of determination R 2 .A value of R 2 = 1 indicates a perfect fit.Fig. 7(b) that it roughly coincides with the centerline of the torus with nonmonotonous behavior as function of κR O and would be still a good approximation, as I have checked, the quality of the fit is a little better if it is used as a free parameter.The main determining factor for describing the full electrostatic potential is, however, ϒ or equivalently .
Finally, in Fig. 7(c), I plot a statistical parameter that quantifies the quality of the surface potential fit, the coefficient of determination or R-squared value R 2 , with R 2 = 1 indicating a perfect fit.From Fig. 7, one can see that the fit is of lesser quality for larger κR O and thicker particles, which is not surprising considering that an infinitely thin torus would be described perfectly by the mapping, and thick particles are a deviation from that.However, the quality is good enough to approximate the full spatial dependence of the electrostatic potential; see also the next section.

APPENDIX C: MORE COMPARISONS OF THE ELECTROSTATIC POTENTIAL FOR AN ELECTROSTATICALLY SCREENED CHARGED TORUS
In Fig. 2 of the main text, I showed some comparisons of the numerically obtained electrostatic potential with the semianalytical approximated ring-charge mapping.Here, I will show that the approximation works even in a relatively "extreme" cases of an even thicker torus R I /R O = 0.8, and extremely thin double layers (up until κR O = 10).The results are shown in Fig. 8. Regarding the charge parameters and shape parameters for these tori, I find for κR O = 1, R = 1.01058, and = 1.1879; for κR O = 5, I find R = 1.0545 and = 1.8078; and finally, for κR O = 10, I find R = 1.0392 and = 2.2751.See also Fig. 7.

APPENDIX D: DERIVATION OF THE TORUS ANISOTROPY FUNCTION WITH YUKAWA MULTIPOLE EXPANSIONS
Sometimes it is beneficial to assess the far-field behavior of the electrostatic potential, which can be obtained with multipole expansion.Such an expansion is useful in deriving, for example, the anisotropy function of a particle; see Eq. ( 8) in the main text.For the dimensionless electrostatic potential of the form a general multipole expansion is given by dr r i 1 ...r i n q(r ).× (r i r j δ kl + r k r l δ i j + r i r k δ jl + r j r l δ ik + r i r l δ jk + r j r k δ il ) with G (0) (r) being the monopolar, G (1) (r) the dipolar, G (2) (r) the quadrupolar, G (3) (r) the octapolar, and G (4) (r) the FIG. 8. Comparison of electrostatic potential of the semianalytical ring-charge mapping with numerical finite-element calculations of the Poisson-Boltzmann equation for a toroidal particle at various values of κR O , along various cuts around the torus as shown in the insets, as well as an enlargement of the comparison between numerical and analytical results close to the particle surface.In all plots, R I /R O = 0.8, which is a thicker torus than the one discussed in the main text (Fig. 1).Furthermore, p / s = 0.2, B /R O = 0.01, and Z = 50.The values of , or equivalently ϒ, and R are determined from Fig. 6.FIG. 9. Comparison of the electrostatic potential of a charged torus as obtained from finite-element calculations, with the Yukawa multipolar-expanded form Eq. (D17) on log-linear scale.The contributions of the various multipoles are highlighted, where monopolar (M), quadrupolar (Q), and hexadecapolar (H) effects are considered.The blue lines are the same as the one discussed in the main text Fig. 1.For the purple lines, the values are cut off because for sufficiently small z at ρ = 0 the potential becomes negative up until quadrupolar order.FIG. 10.Comparison of the numerically obtained anisotropy function with the multipolar expanded anisotropy function obtained from Eq. ( 8).The anisotropy function was obtained for the same parameters as the potential in main text in Figs.1(a), 1(c)-1(e).hexadecapolar basis function tensors, respectively.Note that there are some differences with the unscreened case; for example, the quadrupolar tensor is not traceless, and therefore T (2) cannot be chosen traceless.For line charge distributions, the 2 n -pole moment tensors simplify to In order to compute the tensor contractions for the hexadecapolar term, it is useful to use the formula which is valid because all z components are zero for the ring.Furthermore, the notation G (n) a,b means that x occurs a times and y occurs b times.Specifically, G (4)  i jkl (r)T (4)  i jkl = G (4)  xxxx (r)T (4)  xxxx + 6G (4)  xxyy (r)T (4)  xxyy + G (4)  yyyy (r)T (4)  yyyy , (D15) which results in G (4)  i jkl (r)T (4)  i jkl = Q r R 4   The accuracy of Eq. (D17) is tested against numerics in Fig. 9, showing that at least up until quadrupolar terms are needed to describe the far-field anisotropy sufficiently and hexadecapolar terms increase the accuracy even further at shorter distances from the particle.For r → ∞, where I defined the anisotropy function A(θ ).I compare A(θ ) in Fig. 10 for κR O = 1 for the torus discussed in the main text.
Unfortunately, it is difficult to make comparisons for shorter screening lengths.To appreciate why, consider the example that the anisotropy function from Eq. ( 8) is not sensitive at κR O = 5 for inclusion of the hexadecapolar term; therefore, I can suffice with only the quadrupolar term.In order to numerically compare Eq. (D18) with the asymptotic expression for r → ∞, one must have that 3κR 2 /r 1 (highest order in 1/r quadrupolar term after the constant term).Already at 3κR = 1 the potentials are, however, O(10 −8 ), so one needs to resolve numerically ϕ(r) 10 −8 which is difficult to do considering that the potential close to the surface is O(10 −2 ) for the considered set of parameters; see Fig. 9.
FIG. 3. Semianalytical approximation of the effective pair potential for identical tori [Eq.(10), lines] compared with finite-element calculations (dots) for various particle configurations and screening lengths at B /R O = 0.01 and p / s = 0.2.The colored shaded areas indicate overlap of particle hard cores when the particle distance d is decreased at fixed orientation from d → ∞.Tori have charge Z t = 100 and inner radius R I = 0.5R O .

FIG. 4 .
FIG. 4. Comparisons of the semianalytical approximation of the effective pair potential (lines) with numerical finite-element calculations (dots) between a charged torus and a charged sphere [Eq.(13)] with arbitrarily chosen particle configurations and screening constants.The colored shaded areas indicate forbidden regions due to the particle hard cores when the particle distance d is decreased at fixed orientation from d → ∞.The tori have charge Z t = 100 and inner radius R I = 0.5R O .Spheres have the same surface charge density as the torus with radius a = 0.4R O .For all cases B /R O = 0.01 and p / s = 0.2.

FIG. 5 .
FIG. 5. Comparison of [(a)-(c)] electrostatic potential [Eq.(18)] and (d) the effective pair potential [Eq.(21)] of the semianalytical ellisoidal-shell mapping with numerical finite-element calculations of the Poisson-Boltzmann equation for a triaxial ellipsoidal particle at various values of κa l , with a l being the distance from the origin to the particle surface along the x axis.The potential is compared along various cuts around the ellipsoid, as well as an enlargement of the comparison between numerical and analytical results close to the particle surface.In all plots, p / s = 0.2, B /a l = 0.01, and Z = 100.The shell that has been used to represent the particle has the same ratios of the ellipsoidal main axes, but with a rescalable parameter R l < a l and charge parameter ϒ with the same meaning as for the torus.For precise dimensions of the particle, see the main text.
FIG.7.Result of (a) the charge parameter and (b) shape parameter R as obtained from fits of the surface potential, see Fig.6, for various values of κR O and torus thicknesses as quantified by R I /R O .The results are independent of particle charge Z and Bjerrum length B in the linear screening regime, but depend on p / s (fixed at 0.2 as in Fig.6).In panel (c), we plot a statistical parameter showing the quality of the fits, the so-called coefficient of determination R 2 .A value of R 2 = 1 indicates a perfect fit.
n) (r) is the 2 n -polar basis functions expressed in tensor form and T n is the 2 n -pole moment tensor.Both tensors are symmetric in all indices, and I used the Einstein summation convention.Now consider the DH Green's function G(r, r ) = exp(−κ|r − r |)/|r − r |.In this case, I find up until n = 4,