Rigorous treatment of pairwise and many-body electrostatic interactions among dielectric spheres at the Debye–Hückel level

Abstract Electrostatic interactions among colloidal particles are often described using the venerable (two-particle) Derjaguin–Landau–Verwey–Overbeek (DLVO) approximation and its various modifications. However, until the recent development of a many-body theory exact at the Debye–Hückel level (Yu in Phys Rev E 102:052404, 2020), it was difficult to assess the errors of such approximations and impossible to assess the role of many-body effects. By applying the exact Debye–Hückel level theory, we quantify the errors inherent to DLVO and the additional errors associated with replacing many-particle interactions by the sum of pairwise interactions (even when the latter are calculated exactly). In particular, we show that: (1) the DLVO approximation does not provide sufficient accuracy at shorter distances, especially when there is an asymmetry in charges and/or sizes of interacting dielectric spheres; (2) the pairwise approximation leads to significant errors at shorter distances and at large and moderate Debye lengths and also gets worse with increasing asymmetry in the size of the spheres or magnitude or placement of the charges. We also demonstrate that asymmetric dielectric screening, i.e., the enhanced repulsion between charged dielectric bodies immersed in media with high dielectric constant, is preserved in the presence of free ions in the medium. Graphic abstract


Introduction
The search for methods providing more accurate descriptions and deeper understanding of electrostatic interactions among charged dielectric spheres in ionic solutions remains, despite its long history, a subject of high interest, as proven by a persistent stream of theoretical and experimental publications. This continued interest is motivated by new types of dielectric objects that could be approximated by dielectric spheres and whose properties and interactions are important for other fields of physics, chemistry, materials design and biology.
In order to describe charged dielectric spheres in ionic solutions, it is desirable to choose models for the charged dielectric spheres and the electrolyte solution that may be productively combined, solved consistently, and, ideally, yield a rigorous solution. The most straightforward formulation for the interaction of the charged dielectric spheres is a traditional differential equation boundary value problem of electrostatics. Thus, one has to introduce a compatible model (differential equation) for describing the ionic solution. a e-mail: yyu@ncbi.nlm.nih.gov (corresponding author) It is commonly accepted that, under not too extreme circumstances, each ion in the ionic solution can be considered freely diffusing through the solution in a canonically averaged electric potential φ(r), created by an "atmosphere" of other free ions and ionized colloidal particles. On top of that, the ions are assumed to be thermalized according to Maxwell-Boltzmann statistics. These assumptions give rise to the Poisson-Boltzmann equation: where q s and c s are the charge and average concentration of the ions of species s, o is the dielectric function of the solution, and β is the inverse temperature. It is worth mentioning that the solution of the Poisson-Boltzmann equation, while certainly being useful for providing physical insights, cannot be held as a gold standard of colloidal electrostatics. Indeed, as was originally pointed out by Onsager [1] and Fowler [2], the solutions of the Poisson-Boltzmann equation violate a general reciprocity principle, according to which the work done for bringing an ion j to its position r j in the presence of ion k at position r k must be equal to the work done for bringing the ion k to position r k in the presence of ion j at position r j [3][4][5][6] (see a brief discussion of this in Appendix A). If the average potential φ(r) is sufficiently small throughout the solution, β q s φ(r) 1, the exponential in Eq. (1) can be expanded into a power series and only the first-order term kept, as was originally done by Debye and Hückel [7] (the zeroth-order term vanishes due to the overall charge neutrality): where is the inverse Debye screening length. Using the Debye-Hückel equation (2) instead of the Poisson-Boltzmann equation (1) may in fact be more justified physically as, in the words of Lars Onsager, "... as soon as the higher terms in the Poisson-Boltzmann equation become important, we can no longer expect the ionic atmospheres to be additive, and then the Poisson-Boltzmann equation itself becomes unreliable." [3] Importantly, not only does the solution of the Debye-Hückel equation satisfy the reciprocity principle, it can further be shown [5] that it is the exact small κ limiting form of the Poisson equation for the canonically averaged potential. Thus, curiously, the Debye-Hückel equation has in some ways a firmer physical foundation than the Poisson-Boltzmann equation, in spite of the fact that the former is obtained as a linearization of the latter. Moreover, the range of validity and utility of the Debye-Hückel equation may be expanded by choosing κ appropriately [8][9][10][11]. Since the Debye-Huckel equation obeys certain symmetries violated by the Poisson-Boltzmann equation, allows for rigorous solution, and can have its range of applicability extended by considering κ to be an effective parameter, the Debye-Hückel theory seems a desirable choice for combining with a theory of charged dielectric spheres. There are field theoretical [12] and integral equation [13] techniques for investigating the thermodynamics of ionic solutions that can potentially address the nonlinear regime. However, as we are attempting to combine an electrolyte theory with an electrostatic boundary value problem, and we further anticipate application for simulations, such methods do not seem appropriate for the purposes considered here.
However, one needs to be careful to verify that the Debye-Hückel linearization is valid. In order to do so, one needs to actually solve the electrostatics problem and find the maximum magnitude of the potential, making sure that it is smaller than 1/β = k B T (≈ 26 meV for a solution at 298 K).
In Fig. 1, we plot the electrostatic energy e φ(r) of an ion in an ionic "atmosphere", in the units of k B T , as a function of the distance between the surfaces of two spheres. The potential is obtained as the rigorous solution at the Debye-Hückel level (see Sect. 2) and is calculated at a point where it attains its maximum value, at the surface of one of the spheres, and at the point closest to the other sphere. The parameters of this two-sphere system were chosen to be similar to conditions used in a set of optical tweezer experiments that measured the interaction force between two charged colloidal particles [14]. We see from the figure that for polystyrene particles of radius a = 1 µm, uniformly charged to 10 5 elementary charges e, the Debye-Hückel linearization is valid up to the point where the spheres touch.
The Debye-Hückel linearization is generally less accurate for weaker solutions (smaller screening strength κa), as it is clearly seen in Fig. 1. For demonstration purposes, the smallest κa we used in Fig. 1 was 0.1, while the molarities used in the experiment [14] would result in κa ranging in tens and hundreds. Therefore, the Debye-Hückel equation would still be applicable even if the colloidal particles were even more highly charged.
Another important conclusion that can be drawn from Fig. 1 is that calculating potential at the surface of an isolated sphere in an ionic solution may not be sufficient for judging the validity of the Debye-Hückel linearization. Indeed, the potential is greatly enhanced when two spheres interact with each other. The potential can increase several fold when two charged spheres are in a close proximity, due to significant values of induced surface charge. This circumstance makes finding an accurate solution to the electrostatic problem even more crucial.
Soon after Debye and Hückel proposed their linearized Poisson-Boltzmann model for describing properties of electrolytes [7], the question of how to correctly describe within this model the interaction between charged dielectric bodies was actively debated [15][16][17][18][19].
It was unclear what role the van der Waals interaction plays as compared to the direct electrostatic interaction and even whether electrostatic interaction between similarly charged particles is repulsive or attractive due to the presence of free ions in the solution [20,21].
The consensus model was formulated by Derjaguin and Landau [22,23] and independently by Verwey and Overbeek [24], colloquially joined together and known as the DLVO theory [25]. The DLVO approximation (with or without the ad hoc van der Waals term) is widely used when a quick estimate of interaction energy or force is needed [26]. While the DLVO energy is beautiful in its simplicity and physical clarity and remains in use (see, for a recent example, Ref. [27]), its numerical accuracy and the range of phenomena it can capture are inherently limited.
Beyond the DLVO approximation, there exists a very wide variety of theoretical approaches to colloidal particles in an ionic solution. Not aiming at giving a comprehensive review of these efforts and focusing primarily (but not exclusively) on studies considering dielectric spheres in electrolyte solutions, we can refer the reader to some recent works in the field, employing the method of image charges [28], perturbative meth- Fig. 1 The electrostatic energy e φ(r) of an ion in an ionic solution for interaction of two polystyrene particles of radius a = 1 µm, uniformly charged to 10 5 elementary charges e, calculated near the surface of one of the spheres (r = a) and measured in units of kBT . The ratio eφ(a)/kBT is plotted as a function of separation between the spheres' surfaces for three screening strengths: κa = 0.1 (upper, blue line), κa = 1 (medium, magenta line), κa = 10 (bottom, red line). The dielectric constant of polystyrene is taken to be 2.6, while for the solvent o is taken equal to 78.3. The potential is calculated using the rigorous formalism presented in Sect. 2. Note that condition eφ(a)/kBT 1 is the applicability condition for the Debye-Hückel equation. In the inset, we provide a two-dimensional cross-sectional view of the strength of the potential for the spheres separated by 250 nm, which ranges from approximately 0.11 kBT (red background) to 0.02 kBT (yellow background) ods [29], integral equation methods [30], expansions of induced surface charges into a series of polynomials [31], solving the Poisson equations in specialized coordinate systems [32], discretization of (potentially arbitrary) dielectric surfaces aimed at finding numerical solutions for arbitrary geometries [33,34], direct numerical solution of differential equations [35,36], and so-called hybrid methods [37]. There are, of course, various additional issues such as considering the finite size of the electrolyte ions [38,39] or developing an analytic solution for interaction of spheroidal objects [40], but these are farther afield for the purposes of this article. A number of methods have been developed to describe systems with higher multipole moments. For example, Hoffman et al. [41] considered higher multipoles for a Debye-Hückel system, but only for a single sphere. Also, Boon et al. [42] and de Graaf et al. [43] considered higher multipole moments at the Poisson-Boltzmann level for a single sphere with axially symmetric surface charge. Making progress on the issue of many-body effects, Russ et al. [44] investigated few body forces but were restricted to considering monopole charges and to low order in the number of spheres.
Among the methods using the Debye-Hückel theory and employing an expansion of surface charge into a polynomial series and using the boundary conditions to find the expansion coefficients of such series, we should highlight the rigorous formalism of Fisher et al. [9]. They devised a set of polynomials to serve as the basis set of the expansion and were able to solve the boundary value problem for two identical spheres in an electrolyte solution with either like or opposite point charges at their centers, thereby obtaining an expression for the interaction energy. Even though a direct analytic comparison would be cumbersome, their result for this symmetric two-sphere system is numerically equivalent to that obtained within the general formalism developed recently [39] (even the convergence rate is the same, see Fig. 9 later).
Sushkin and Phillies [45] attempted to address the general case of two arbitrary spheres with an arbitrary collection of point charges within each. Their formalism is physically sound, but suffers from unpolished presentation, calling for a great deal of further simplification and clarification. Moreover, the printed manuscript probably contains typos. 1 During the recent decade, Besley and Stace et al. have achieved significant progress in various aspects of inter-actions between two dielectric bodies in vacuum and in an electrolyte medium, reported in a series of papers [40,[46][47][48][49][50][51]. Their formalism [46], although imposing azimuthal symmetry with respect to the axis connecting the centers of the two spheres, was demonstrated to provide useful insights, such as dependence of attraction between two spheres with same charge immersed in a medium of lower dielectric constant on asymmetry of charge and radii of the interacting spheres [47,50]. Their analysis uses a re-expansion of modified Bessel functions about a new center to allow matching of the boundary conditions at the surface of each sphere. However, multiple re-expansions are required in order to obtain the necessary Legendre polynomials, making controlling the numerical accuracy somewhat difficult since rather than a single maximum multipole max , the multiple re-expansions introduce additional sums that need to be cut off in a consistent manner.
In a series of efforts [52][53][54][55], we have undertaken the development of a rigorous method for determining to any desirable precision the interaction energies for an arbitrary number of dielectric spheres in a dielectric medium. The most recent work [39] contains a general formalism, rigorous within the Debye-Hückel approximation, for describing interactions among dielectric spheres immersed in an electrolyte solution. The only other surface charge method of a similar generality and rigor that we are aware of is due to Lotan and Head-Gordon [56]. The downside of their method, however, is that transformations of coordinates when reexpanding at different locations are handled via an iterative numerical procedure. This significantly complicates computations, reduces their efficiency, and, even more importantly, introduces another source of uncertainty in the accuracy of the numerical results on top of cutting off polynomial expansions of the surface charge distributions.
In this article, we apply the general formalism [39] to detailed analysis of interactions in the case of two, three and four spheres. A terse summary of the general formalism may be found in Appendix B. We work out the simplifications that are possible for the particular case of two spheres and even further simplifications for the case of two spheres with axially symmetric charge distributions. In the latter case, it is possible to derive a simple expression for the force in terms of the same set of linear equations as for the energy. The significant new physical results are as follows. We demonstrate in Sect. 2.1 that asymmetric dielectric screening [52,57], i.e., the enhanced repulsion between charged dielectric bodies immersed in media with high dielectric constant, known to occur in the absence of ions, is preserved in the presence of free ions in the medium. In Sect. 3, we compare the baseline DLVO approximation [25] to the results of the rigorous formalism. Indeed, only by comparing to the exact solution of the Debye-Hückel problem can the accuracy of DLVO theory (or any of its proposed modifications or replacements) be correctly and reliably assessed. We highlight the circumstances under which one can expect larger errors for the DLVO theory. In the last section, we discuss the accuracy of approximating the full interaction energy as a sum of pairwise interactions using examples of three-and foursphere systems. Finally, we illustrate how the speed of convergence depends on variation of radii, charge magnitude and charge placement in the system.

Two spheres
In order to represent charged colloidal particles in an ionic solution, we consider a system of dielectric spheres in a medium (water) with dielectric constant o containing freely diffusing, thermalized ions. Sphere i has dielectric constant i , radius a i , position r i , and free charge that can be represented by standard electric multipole moments q lm (l = 0, 1, . . . , ∞ and m = −l, −l + 1, . . . , l) or equivalently [39,54] by surface moments Q lm . This free charge distribution is fixed (we do not consider a dynamic chemical charge regulation scheme for the free charge), but there will be additional induced charge at the interface between the regions of differing dielectric constant. The spheres are embedded in an electrolyte solution that will be modeled by the Debye-Hückel theory and whose properties are therefore summarized by the inverse Debye screening length κ. We will, for the most part, regard κ as a parameter. Thus, Eq. (3) will not be binding.
Within the framework of the general formalism for this system developed in Ref. [39] and briefly recapitulated in Appendix B, the distribution of the fixed (also called free) charge ρ(s) at each sphere is expanded into a series using the basis set of spherical harmonics Y m .
The (yet unknown) induced charge distribution at each sphere is also expanded with the same basis set, and the boundary conditions on the potential at the surface of each sphere then lead to a system of linear equations. The variables in this system are scaled spherical components of the net (free + induced) charge distri-butionQ + m . Having determined these components by solving the system of equations, one can subsequently find the induced charge distributions on each sphere, the electric potential at any point of space, as well as the electrostatic interaction energy.

Two spheres with arbitrary charge distributions
For a system of only two spheres, the general system of linear equations simplifies not only due to fewer expansion centers, but also due to the fact that for two spheres there is only one vector connecting the centers of the spheres, and this vector L ≡ L 1→2 = −L 2→1 can always be directed along the z-axis. The spherical harmonics depending on the orientation of this vector are, therefore, only nonzero for zero momentum projections: As small as it seems, this circumstance allows one to derive a quite simple, yet still rigorous, system of equations for solving the electrostatics problem of two spheres in an ionic solution.
With (4) in mind, the general system of equations (39) can be reduced to: Here for convenience we use explicit enumeration of the spheres. The source terms in (5) are the spherical components of the free charge distribution (surface moments) Q m , where a is the sphere's radius and q m are regular multipole moments of the charge distribution ρ(s).
As is expected for a spherical harmonics expansion for a Helmholtz equation, the radial coefficients are the modified spherical Bessel functions: i (z) are the modified spherical Bessel functions of the first kind with i (z) = (i) − j (i z), while k (z) are the modified spherical Bessel functions of the second kind with The functions are calculated at the dimensionless radius of one of the spheres κa 1 or κa 2 .
The functions K( , x) and I( , x) are introduced solely for the brevity of notation; they depend on the dielectric constant of the medium o and on the dielectric constant of the corresponding sphere: Finally, the function H m 1 (κL) is defined as: it is a reduced version of a more general factor H m 1m1 (κL k→j ), Eq. (36), that in the general formalism takes care of the mutual location and orientation of different spheres. Note, the quantity H m 1 (κL) is symmetric with respect to the interchange ↔ 1 , allowing for a more efficient numerical implementation and further simplifications.
With the above definitions, the system of linear equations (5) with an appropriately chosen maximum number of components max can be solved, and the unknown components of the scaled net charge distributionsQ 1+ m andQ 2+ m can thus be determined. Now, with the knownQ 1+ m andQ 2+ m , the interaction energy can be found as: Note that Eqs. (5) and (9)  m can also be used to determine the electrostatic potential between the spheres as: Here r 1 is the point where the potential is to be found. We chose the system of coordinates to be associated with the first sphere, but the second sphere's coordinates can also be used with the appropriate change of indices. The direction of the vector r 1 can be arbitrary, but the distance where Eq. (10) is valid is limited by the inequality a 1 < r 1 < L − a 2 .
As an illustration of capabilities of this general method, let us consider electrostatic interaction in a non-axially symmetric system. In Fig. 2, we plot interaction energies between two spheres with physical dipoles oriented perpendicular to L 1→2 . The dipoles are created by two point charges, shifted from the corresponding sphere's center by half of its radius. The charges are equal in magnitude and opposite in sign, so the net charge of each sphere is zero. The interaction energy is calculated for three different orientations of these dipoles: parallel, anti-parallel and orthogonal. As expected, the spheres with parallel dipoles repel and the spheres with anti-parallel dipoles attract. There are two noteworthy features, however, in Fig. 2. The first one is the asymmetry in repulsion of parallel dipoles and attraction of anti-parallel ones. This stems from asymmetric dielectric screening [52], a known effect in the interaction between dielectric objects whereby either repulsion (if in < out ) or attraction (if in > out ) is enhanced as compared to interaction of point charges. For a discussion and an application of this effect to interaction of dielectric spheres, see, e.g., Refs. [50,57]. The second noteworthy feature is the fact that the interaction between two orthogonal dipoles (shown in magenta) is not zero, as one might have anticipated. Indeed, if the dielectric spheres were absent, then the interaction would vanish. However, the asymmetric dielectric screening is an independent feature due to the dielectric mismatch between the spheres and the solvent and persists even when the interactions of the point charges cancel.

Two spheres with axially symmetric charge distributions
When axial symmetry is lacking, such as in the example of two dipoles considered above, all the components of the scaled free and net charge distributions, Q m and Q + m , can be nonzero. This means that one has to solve the complete system of 2( max + 1) 2 linear equations (5). However, in the case of axial symmetry, only the m = 0 components Q 0 andQ + 0 will be nonzero. The new system will only have 2( max + 1) equations, and the calculations will simplify drastically.
System (5) is then reduced to: For further simplification, one can use explicit expressions for Clebsch-Gordan coefficients with zero momentum projections [58] to obtain H 0 1 (κL) in terms of factorials, where 2g = + 1 + 2 must be even. Moreover, it is now easy to find the force acting on each sphere, because for charge distributions axially symmetric around the vector L, the force is directed along L. In this case, in order to determine the force, we can take a derivative of the interaction energy with respect to the distance L (instead of taking a gradient as in the most general case): For computational purposes, one can again express the derivative of the modified spherical Bessel as in Eq. (41): The derivatives of the components of the scaled net charge distributionsQ + 0 can be found by taking derivatives of both parts of Eq. (11): Assuming that the system of linear equations (11) has already been solved and introducing quantities we obtain a system of linear equations for the negative derivatives: Note that systems (11) and (17) have the same matrix of coefficients; only the source terms are different. This makes a computational implementation of force calculation much more efficient. This increased computational efficiency, appearing in the case of axially symmetric systems due to the reduced size of matrix (11) and the possibility of its re-purposing for determining the force, allows calculating electrostatic interactions in a rigorous fashion down to very short distances, almost to the touching point, where the number of contributing terms may increase to hundreds.
Janus particles (see, e.g., a recent paper [59]) could serve as a practical example of such an axially symmetric system. Within the present formalism, they are easily modeled as two patches of positive and negative charge, uniformly distributed over a spherical cap of a polar angle θ 0 . In this case, the spherical components of the free charge are given by: where P (x) are Legendre polynomials. It is interesting to note that for reasonable values of the opening polar angle θ 0 the scaled free charge components are close to those of a simple dipole oriented along the z-axis, as can be seen from the small-θ 0 expansion (19) Our numerical results confirm this conclusion.
Another instance of a directly observable situation where an axially symmetric charge distribution might prove useful is modeling interaction of two dielectric beads that are attached to a support. For example, in an optical tweezers experiment one of the colloid particles is held in place by a pipette [14], which makes an otherwise spherically uniform charge distribution only axially symmetric. (Our calculations show that the presence of a pipette in these kinds of experiment would be important for low κa regimes.)

Two spheres with point charges at their centers
Finally, let us write down the expressions for the interaction energy and the force for the important case when both spheres have only one point charge located at their center. Let's denote these charges q 1 and q 2 . Then, according to Eq. (6), The interaction energy now has two simple terms centered at each sphere and only one summation in the cross term. Indeed, only the = 0, m = 0 components in the first two sums in Eq. (9) survive, while the compounded sum over , m and 1 in the third term turns into two independent sums over and 1 (which can be combined by renaming 1 back to ).

Comparison with the DLVO theory
The presented formalism is exact and therefore allows a quantitative characterization of errors arising from using the DLVO approximation. Within the DLVO approximation, the electrostatic interaction energy is given by the term: where D = L − a 1 − a 2 is the separation between the spheres' surfaces. Generally, the accuracy that can be achieved with the DLVO theory is governed by the number of terms in the polynomial expansion of the surface charge that contribute significantly to the interaction energy. For large separations and/or very short screening lengths, the = 0 term dominates and the DLVO approximation is accurate. For short separations and moderate screening lengths, however, the DLVO approximation can lead to large errors. In Fig. 3, we plot the ratio of U DLVO int to the full interaction energy U int given by Eq. (21) for three dimensionless screening strengths κa = 0, 0.1, 1. As the strength of screening increases, the DLVO approximation becomes worse at small separations of the spheres, while getting better at large separations. The accuracy of the DLVO theory also suffers for asymmetric charge distribution and asymmetric radii (compare the solid and the dashed lines in the figure). Indeed, for two identical spheres only the first two terms in the expansion are enough to achieve accuracy of one per cent even at the point where the spheres touch, the configuration that presents the most stringent test of accuracy and convergence. However, for highly asymmetric radii and charge distributions one may need a hundred terms to get to the relative error of 10 −3 .
While the inaccuracies of the DLVO theory in interaction energies are usually moderate, they are much higher for the interaction forces, the physical quantities that are directly measured by atomic force microscopy or optical tweezers experiment, see, e.g., Refs. [14,60]. As an illustration, in Fig. 4 we plot the ratio of the DLVO force to the exact, fully converged interaction force of Eq. (22), for the two systems for which the ratio   . 4 The ratio of the approximate DLVO interaction force to the fully converged force of Eq. (22) in the plane of surface-to-surface separation D and screening strength κa for the two identical spheres (inset) and for two spheres with different radii and charges. The parameters of the system are the same as in Fig. 3. The contours represent (going from the top right corner toward the bottom left corner) the ratios of 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, and 0.2 of energies is plotted in Fig. 3. For a more comprehensive perspective, we plot the force ratio in the plane of separations D and screening strengths κa. The region of the parameter space for which the ratio of the DLVO force to the exact force is greater than 0.9 (DLVO accurate within roughly 10% error) has a white background. The dark red background corresponds to fivefold errors of the DLVO approximation (the 0.2 ratio). For large separations D and values of κa (the top right corner), the ratio of the DLVO's result to the exact one is nearly one, while the ratio decreases (for the asymmetric system-dramatically), toward the area of small separations and weak screening.
We also present interaction forces for conditions similar to those of experiments measuring forces between two colloidal particles [14], as shown in Fig. 5. For a 1 µm particle charged to 10 5 elementary charges (0.008 elementary charges per nm squared), at small separations the forces again differ by several fold, while the interaction energies differ by only several tens of per cent.
We conclude, therefore, that the DLVO approximation provides accurate interaction energies and forces for large separations and for strong screening, while for short separation, weak screening and for asymmetric charge/radii combinations the DLVO theory predictions can be misleading. This is also true for systems involving multiple interacting spheres, where manybody effects, considered in the next section, can play an important role.

Many-body effects
Obtaining a general solution to the problem of multiple dielectric spheres in an electrolyte, not restricted by any symmetry considerations, is a much harder problem, as compared to calculating the interaction between just two spheres. Not only does the general solution require more computational resources, its results are more difficult to predict and analyze. That is why it is so tempting to approximate the interaction within a set of spheres by the sum of their pairwise interactions. The difference between the results obtained within the pairwise approximation and within a full formalism is frequently referred to as "many-body effects". In any particular physical model in which many-body effects are possible, the role they play might range from vital to, on the contrary, a negligible one. In this section, we use the general formalism [39] to study the importance of many-body interactions for systems of three and four dielectric spheres in an electrolyte.
Earlier it was demonstrated [61] that many-body effects are important in interactions of several dielectric spheres in vacuum, and their significance increases with the number of interacting spheres. For example, if three positively charged spheres were arranged in an equilateral triangle, the charge density on each sphere had a broad distribution of negative surface charge centered in the direction of the triangle's center. On top of that broad distribution, there were two additional smaller peaks oriented toward the other two spheres. However, within the pairwise approximation the negative surface charge distribution was less pronounced and, moreover, had a simple singly-peaked sinusoidal shape, a result of a superposition of two sinusoidal pairwise charge distributions. Thus, the pairwise approximation in this example failed both in predicting correct magnitude of the surface charge density and in predicting possible fine features of the surface charge distribution. Such inaccuracies in the surface charge density led to substantial errors in the interaction energies found within the pairwise approximation [61]. Let us now examine how the presence of mobile ions in the dielectric medium influences the significance of the many-body effects.
From the physical viewpoint, it is clear that at large distances, κL 1, pairwise interactions should provide a very accurate estimate of the total interaction energy and the shorter the Debye length becomes, the more accurate this estimate gets. This is due to the fact that at large distances electrostatic interactions are dominated by the = 0 term, while the contribution of higher-order terms is small and, moreover, gets screened more and more effectively with shorter Debye length. With respect to another length scale, κa, in the limit of a very short Debye length, κa 1, the direct Coulomb interaction term = 0 would also dominate even as the spheres touch, since the mobile ions would completely screen the presence of the other spheres, suppressing the manifestation of many-body effects. However, it is not obvious what happens in a situation when the Debye length is comparable to the radii of the spheres, κa ∼ 1, and the spheres are not too far apart, a ∼ L.
Our calculations show that in this regime, the presence of mobile ions in the medium actually enhances Fig. 6 The ratio of the interaction energies between three spheres forming an equilateral triangle with the side L, calculated as a solution of the consistent many-body Eq. (44) and as a sum of three pairwise interactions between two spheres. The ratio is calculated for point charges at the center of each sphere (black lines marked by a gray geometry sketch) and for point charges shifted by half the radius from the center of each sphere toward the geometric center of the system (green lines and the green geometry sketch). Three dimensionless screening strengths κa = 0.1 (thick lines), κa = 1 (thin lines) and κa = 2 (dashed thin lines) were used. The charge q and the radius a are one in atomic units; the dielectric constants inside the spheres are 1 = 2 = 4, and o = 80 for the medium. The cutoff values of max = 20 (shifted charges) and max = 15 (centered charges) were used to achieve convergence to the relative accuracy level of 10 −3 the importance of collective effects, many more higher multipoles contribute significantly to the interaction energy, and the quality of the pairwise approximation gets much worse.
This phenomenon is clearly seen in Fig. 6. In this figure, we plot, for a system of three spheres, the ratio of interaction energies found within the pairwise approximation (the sum of three identical pairwise interactions between the spheres) and within the rigorous formalism, given by Eqs. (39) and (44). The spheres form an equilateral triangle. Inside each sphere, there is a positive point charge that we position either at the sphere's center or shift by a half of the radius toward the geometric center of the triangle. The ratio is plotted for three values of the screening strength, κa = 0.1, 1, 2.
As anticipated, at large separations and as κ increases, the curves move closer to 1, i.e., the pairwise approximation works better. However, at shorter distances the dependence on κa is not trivial. Consider the geometry where the spheres touch (L = 2a) and the charges are at the centers of the spheres. For κa → 0, the error of the pairwise approximation is slightly more than 4%. It increases with κa, reaching about 10% when κa ∼ 1 and then gradually decreases again. 2 These errors are much bigger for the case when the point charges are shifted off the spheres' centers toward the triangle's center. The reason is that in this case the point charges are closer to the surfaces where the spheres touch, so the surface charge distributions are more peaked and more higher-order terms become important. In contrast, if the point charges are shifted away from the triangle's center, the surface charge distribution is peaked at the "back side" of each sphere, with a wide, featureless distributions on the "front side" of the spheres, where they touch. In this case, the pairwise approximation becomes more accurate (not shown in Fig. 6). This observation is in line with the conclusion made for interactions of several spheres in nonionic solutions [61], where the direction and amplitude of peaks in the surface charge density affected the accuracy of the pairwise approximation.
We observed the same behavior of the many-body effects also for a system of four spheres, see Figs. 7 and 8. The spheres are arranged in a tetrahedron, so the distances between their centers are equal. Each sphere can have either a point charge at its center or two point charges, equal in magnitude, but of opposite signs. The latter dipoles are oriented in two spheres toward the geometric center of the tetrahedron and away from it in the other two spheres.
As in the case of three spheres, for large distances the pairwise approximation is more accurate at larger κ, while at shorter distances the many-body effects are important and the pairwise approximation can greatly Fig. 7 The interaction energies in atomic units between four spheres forming a tetrahedron with the side L. Each sphere (of a radius a) has a point charge q at its center. The interaction energies are calculated as a fully converged solution of Eq. (44) (red lines, marked by a tetrahedron icon) and as a sum of six pairwise interactions between two spheres (blue lines marked by a two-sphere icon). Three dimensionless screening strengths, κa = 0.1 (thick lines), κa = 1 (thin lines) and κa = 2 (dashed thin lines), were used. The charge q and the radius a are one in atomic units, the dielectric constants inside the spheres are 1 = 2 = 4, and o = 80 for the medium Fig. 8 The ratio of the interaction energies between four spheres forming a tetrahedron with the side L, calculated as a fully converged solution of Eq. (44) and as a sum of six pairwise interactions between two spheres. The ratio is calculated for point charges at the center of each sphere (black lines marked by a gray sphere) and for physical dipoles in each sphere (green lines marked by a green sphere). The physical dipoles are formed by two point charges, equal in magnitude and opposite in sign, shifted from the cen-ter of the sphere by half radius, so the total length of each dipole is a, the radius of the sphere. Two of the dipoles are pointed toward the geometric center of the system by their positive end, while the other two dipoles-by their negative end. Three dimensionless screening strengths, κa = 0.1 (thick lines), κa = 1 (thin lines) and κa = 2 (dashed thin lines), were used. The charge q and the radius a are one in atomic units, the dielectric constants inside the spheres are 1 = 2 = 4, and o = 80 for the medium Fig. 9 Convergence of interaction energy for two identical spheres of radii a and charge q located at the touching point L = 2a (squares) and with a gap of one radius between them, L = 3a, (circles). Open squares and circles are obtained within the present formalism, while filled circles and squares are obtained using Fisher's formalism [9].
The screening strength is κa = 0.1, the charge q and the radius a are one in atomic units, the dielectric constants inside the spheres are 1 = 2 = 4, and o = 80 for the medium. The relative errors are calculated with respect to the fully converged value obtained for max = 100 Finally, we note that the correct account of manybody effects can only be achieved if enough terms are included in expansion of the potential/surface charge. Unfortunately, the convergence can be quite slow when spheres' surfaces approach each other. This general characteristic of methods that rely on re-expanding spherical Bessel functions and Legendre polynomials at shifted origins is illustrated in Fig. 9. There we plot convergence of interaction energy for two identical spheres at the touching point and with a gap of one radius between the spheres. The convergence is plotted both for our method and for Fisher's method [9], which employs a different set of basis polynomials. However, the rate of convergence is nearly identical for this system of identical spheres.
Apart from the dependence on the distance between the spheres, the convergence is generally slower for systems with large variations in the relative magnitudes of radii and charges (see, e.g., [50]) or with large variations between charge distributions within each sphere (for instance, when a point charge is located not in the center of a sphere, but close to its surface).
In Figs. 10 and 11, we illustrate this point by plotting the relative error of the interaction energy as a function of cutoff value max for the following pairs of systems: two identical spheres versus two spheres with different radii and charges (the system from Fig. 3), three spheres with point charges at the center versus three spheres with point charges off center (the system from Fig. 6), four spheres with point charges at the center versus four spheres with zero net charge, but with physical dipoles inside (the system from Fig. 8).
For each of these pairs, we see that variation in the charges, radii, or off-center placement of charges significantly reduces the rate of convergence. In some cases, the rate of convergence fluctuates even for large max . For four spheres, these irregularities span only one or two expansion orders, while for three spheres we observe a more regular pattern in the fluctuation of convergence rate. We have carefully examined possible sources of numerical errors and, using arithmetic of arbitrary precision and various methods for solving the system of linear equations, verified that these peculiarities seem to be a result of a complex interplay between contributions of many different terms of the expansion, each of which has its own convergence rate and relative importance as a function of max . Further study of the convergence patterns is needed, but care should be taken when calculating interaction energies for asymmetric systems.

Appendix A: Violation of reciprocity in the Poisson-Boltzmann framework
To provide context, we briefly review the discussions by Onsager [3] and McQuarrie [6] of the Poisson-Boltzmann and Debye-Hückel equations. For a neutral solution of positive and negative ions in water ( o), these differential equations provide < φ(r, r1) > (1) , the canonically averaged electric potential at r while keeping a single ion fixed at r1.
Taking the square gradient of the expression for the canonical average of φ gives the first equality. The second equality is by the definition of the radial distribution function g1s(r, r1), where qs and cs are the charge and average concentration of the ions of species s and β = 1/kBT is the inverse temperature. The final step is the approximation that gives the well-known Poisson-Boltzmann equation: where we have without loss of generality set r1 = 0 and used streamlined notation for the potential in accordance with common practice. Expanding the exponential to first order and noting that the zeroth-order term vanishes due to charge neutrality gives the Debye-Hückel equation: where κ 2 = 4π β o s csq 2 s . However, we are not wedded to this expression for κ. As noted by Fisher et al. [9], the range of validity and utility of the Debye-Hückel equation may be expanded by considering κ to be a system parameter. More recently, Janacek and Netz [10] have performed numerical studies that suggest the use of an effective screening parameter κ eff . The work of Alexander et al. [8] also suggests the expanded applicability of Debye-Hückel theory when using an effective κ.
Let us now consider some symmetry requirements. Generalizing the above notation so that < φ(r k , rj) > (j) is the canonically averaged electric potential at a point r k while keeping an ion fixed at point rj (and vice versa). By imposing translational symmetry and particle permutation symmetry, the radial distribution function satisfies g jk (rj, r k ) = g jk (|rj − r k |) = g kj (r k , rj). Note that the potential of mean force w is defined via e −βw j,k (r j ,r k ) ≡ g jk (rj, r k ), implying w j,k (rj, r k ) = w k,j (r k , rj) .

(29)
The reciprocity principle requires that the energy of the system be independent of the order in which it is assembled. Thus, one may bring ion k into position with ion j already fixed or vice versa. The left-hand side above is the work required to bring ion j into position with ion k already fixed, while the right-hand side corresponds to work required for bringing in ion k with ion j already fixed. Expressing this in terms of the separation and taking the negative gradient, one sees that this is also an average form of Newton's third law.
The approximation from (25) to (26) thus implies the need of having qj < φ(rj, r k ) > (k) = q k < φ(r k , rj) > (j) , (30) which is unlikely to hold when nonlinear terms are included. As Onsager stated, we have to expect the left-hand side differs from the right-hand side "except in very symmetric cases." [3] Interestingly, it is straightforward to verify that the solution the the Debye-Hückel equation satisfies this requirement. Moreover, this reciprocity has been explicitly proven [39] for the application of Debye-Hückel theory to an arbitrary collection of dielectric spheres in solution.
Another symmetry requirement arises in the form of a reciprocal relation. Since U = (1/2) j qj i qi/( o|ri − rj|), it is easy to see that:

∂F ∂qj
= φj (31) where F is the Helmholtz free energy and φj is the average potential at rj (excluding the contribution from qj). Because the total differential of the electrostatic portion of the Helmholtz free energy is given by: equality of the mixed partial derivatives implies: Importantly, the solution to the Debye-Hückel equation satisfies these requirements, but solutions to the Poisson-Boltzmann equation will not. Furthermore, it can be shown [5] that the Debye-Hückel equation is the exact small κ limiting form of Eq. (24). Evidently, the Debye-Hückel equation has in some ways a firmer physical foundation than the Poisson-Boltzmann equation in spite of the fact that the former is obtained as a linearization of the latter. Moreover, as Onsager [3] points out, "as soon as the higher terms in the Poisson-Boltzmann equation become important, we can no longer expect the ionic atmospheres to be additive, and then the Poisson-Boltzmann equation itself becomes unreliable." A brief comment on van der Waals forces is perhaps appropriate at this point. The van der Waals interaction arises at the atomic scale as a result of quantum mechanical fluctuations. Thus, the ions in an ionic solution might well experience van der Waals interactions. However, while these forces are significant for atomic scale objects such as the ions in the solution, it will generally not be the case for the much larger colloidal objects (macroions), which deserve a classical approach. Nevertheless, this is not to say that there are no interactions of the form r −6 among colloidal particles. An expansion in inverse powers of the separation r of the electrostatic interaction of polarizable objects in a polarizable medium will certainly include terms such as r −6 [61], but their origin lies in continuum macroscopic electrostatics.

Thus,
Finally, let us note that in the limit of κ → 0 (no mobile ions in the solvent), the presented formalism reduces to the formalism of dielectric spheres in a dielectric medium [55], as shown in Ref. [39].