Exact polarization energy for clusters of contacting dielectrics

Resolving the surface charges is needed to evaluating the electrostatic interactions among particles in close proximity, and is challenging because polarization is intrinsically a many-body effect, depending on the positions of all particles. For instance, careful measurements in colloidal suspensions have shown that the inter-particle force is non-additive, which is at least partially due to the polarization effect.⁴ In another set of experiments on metallic nanoparticles, it has been found that the aggregation of multiple particles surrounding a charged particle can be stabilized by the polarization effect alone.^{3, 5} More dramatic demonstration of such polarization effects is found in the so-called like-charge attraction caused by the strong polarization effect, for particles of large size or permittivity ratios.^{6, 7}

Computational methods, such as the boundary element method,⁸ the spectral methods ^{9, 10} and the image method,¹¹ have been developed to account for this polarization effect. In practice, these methods all need to evaluate the induced surface charge densities in one form or another, and have been successfully applied to study a wide range of problems involving the aggregation of polarizable particles.^{12, 13, 14} However, when particles are in close proximity, the surface charge densities apparently diverge because the electric field in the narrow gap region is strong even for a small difference in the electrostatic potentials of particles. This is analogous to the divergent lubrication force between approaching solid particles in the Stokesian regime.¹⁵ Resolving these apparently diverging charge densities requires a high spatial resolution that becomes computationally prohibitive for nearly touching particles.⁸ For conductors like metalic particles, the functional form of the divergent charge density has been identified and used to isolate the singularity obscuring numerical calculations, and to obtain the exact energy and force for contacting particles.⁵ For dielectric particles, the question remains unsolved.

One way to reveal these singular surface charges is to consider two conducting particles separated by a small gap, as sketched in Fig. 1. The points on the two surfaces separated by the minimum distance, or gap distance $h$, are referred to as contact points. The surface charge density is proportional to the strength of the normal component of electric field on surfaces, according to Gauss’s law. When $h$ is small, the electric field lines in the gap region are nearly parallel to the line connecting the two contact points. The strength of the electrical field is, in turn, proportional to the difference in the surface electrical potentials. This relation between the surface potential and the electrical field is analogous to that between the longitudinal velocities of converging particles and the transverse velocity of squzzed fluid velocity in the lubrication theory.¹⁵

To calculate the surface charge density, both the surface potentials and the vertical separations are needed. Since conductors are equipotential, the surface potentials can be set to $V$ and $0$ respectively. The vertical separation depends on the radial distance $r$, and can be expressed within the Derjaguin approximation as $h+r^{2}/a$, where $a$ is the harmonic average of radii of curvature $a_{1}$ and $a_{2}$ at the two contact points, i.e., $a^{-1}=(a_{1}^{-1}+a_{2}^{-1})/2$. Here, the spherical apexes are assumed for simplicity, but more general curvatures can be treated similarly. Consequently, the field strength at radius $r$ is $E(r)=V/(h+r^{2}/a)$, and the surface charge density $\sigma(r)=\frac{1}{4\pi}\epsilon_{\rm out}V/(h+r^{2}/a)$ on the top surface, where $\epsilon_{\rm out}$ is the medium permittivity. Integrating $\sigma(r)$ for $r$ from $0$ to $R_{0}$, where $R_{0}$ is a cutoff of order $a$, results in the singular part of the surface charge

$\displaystyle Q_{\rm s}$

$\displaystyle=\frac{V\epsilon_{\rm out}}{4\pi}\int_{0}^{R_{0}}\!{\rm d}r\frac{2\pi r}{h+r^{2}/a}\simeq\frac{V\epsilon_{\rm out}a}{4}\ln\left(\frac{a}{h}\right).$

(1)

The unity in the logarithmic term is dropped because $R_{0}\simeq a\gg h$. Further, the difference between $R_{0}$ and $a$ is neglected because it only leads to a constant shift. Equation (1) shows how surface charges of the upper surface become singular as the gap distance decreases. That of the lower surface is also singular, but of opposite sign. In the limit $h\to 0$, the energy does not blow up because the potential difference $V$ vanishes once particles form contact.

The ratio of the singular charge $Q_{\rm s}$ and the potential difference $V$ gives the singular part of the capacitance $c_{\rm s}=\epsilon_{\rm out}a\ln(a/h)/4$, which has been found previous for spherical dimers.^{6, 7} The above analysis shows that this singular capacitance is local. Thus, the same singularity applies to non-spherical particles, and for aggregates of multiple particles. This fact has been employed to find the exact energy for clusters of contacting, conducting particles.⁵ In this work, the full capacitance array is first numerically calculated for an ensemble of conducting particles at finite but small separations. The singular contribution is then subtracted, leaving a regular part that can be extrapolated to $h=0$. Finally, when the singular and regular parts are pieced together, the variation of energy with separation is obtained.

For the dielectric case, a straightforward generalization of the above treatment fails, because the particles are not equipotential. The potential difference $\Delta V(r)$ needed to evaluate the charge density in eq. (1) can not be fully specified by the potential difference between the contact points $\Delta V(0)$. The variation of $\Delta V(r)$ with $r$ in the gap region is expected to be quadratic, but the curvature is unknown a priori. In an early, and analogous work on thermal conduction of composite materials, Batchelor et al.¹⁶ noticed that the potential distribution, which determines the surface charges, is itself dominated by the contribution from the surface charges nearby, so that a self-consistent treatment is needed.

To illustrate this point, we consider the dielectric case in Fig. 1. Let the surface potentials of the two particles be $\phi_{1}(r)$ and $\phi_{2}(r)$, and the surface charges be $\sigma_{1}(r)$ and $\sigma_{2}(r)$. The surface potential $\phi_{i}(r)$, with $i=1,2$, can be calculated by integrating the coulombic potential of the respective surface charges. Near the contact point, we have

$$\phi_{i}(r)=\tilde{V_{i}}+\frac{1}{2\pi\epsilon_{i}}\int_{0}^{\infty}\!\!{\rm d}r^{\prime}\!\!\int_{0}^{2\pi}\!\!{\rm d}\theta^{\prime}\frac{r^{\prime}\sigma_{i}(r^{\prime})}{\sqrt{r^{2}+r^{\prime 2}-2rr^{\prime}\cos\theta^{\prime}}}.$$

(2)

Here, $\tilde{V}_{i}$ are the contributions to the surface potential from the surface charge outside the contact region. The integral are those from the surface charges in the contact region. The prefactor is $\frac{1}{2\pi}$ (instead of $\frac{1}{4\pi}$) because of the well-known jump condition for surface potentials.¹⁷ The distance at denominator is approximated using that for the flat surface, which leads to negligible error because only the contact region is of concern. The upper bound is set to infinity for convenience; as we shall see below, the singular surface charge density dies off rapidly outside the contact region.

The surface charge density $\sigma_{i}(r)$ in eq. (2) are related to the potential difference, $\Delta V(r)\equiv\phi_{1}(r)-\phi_{2}(r)$. By analogy to the conductor case, the electric field is approximately vertical and its magnitude is $\Delta V(r)/(h+r^{2}/a)$. Then applying Gauss’s law, we get the charge density on the top surface

$$\sigma_{1}(r)=\epsilon_{\rm out}\left(1-\epsilon_{{\rm r},1}^{-1}\right)\frac{\Delta V(r)}{h+r^{2}/a},$$

(3)

where $\epsilon_{{\rm r},1}\equiv\epsilon_{1}/\epsilon_{\rm out}$. The charge density $\sigma_{2}$ is given similarly, but with a negative sign. Substituting the charge densities to eq. (2) and taking the difference gives an integral equation for $\Delta V(r)$. The dependence on the unknown, $\tilde{V}_{1}-\tilde{V}_{2}$, can be factored out by introducing an auxiliary, $f(r)\equiv 1-\Delta V(r)/(\tilde{V}_{1}-\tilde{V}_{2})$, which satisfies

$$f(r)=\frac{\epsilon_{{\rm r},1}^{-1}+\epsilon_{{\rm r},2}^{-1}}{2\pi}\int_{0}^{\infty}\!\!{\rm d}r^{\prime}\ \frac{1-f(r^{\prime})}{h+r^{\prime 2}/a}\frac{4r^{\prime}}{r+r^{\prime}}K(x).$$

(4)

Here the integral over the azimuthal angle has been replaced with the complete elliptic function of the first kind $K(x)$, where $x\equiv\frac{4rr^{\prime}}{(r+r^{\prime})^{2}}$. When $\epsilon_{1}=\epsilon_{2}$, eq. (4) reduces to Batchelor’s original result, eq. (4.5) in ref.16. Since $f(r)$ is proportional to the contribution from the surface charges in the contact region, we see that it is dominated by particles of lower permittivity. When both permittivities approach infinity, we have $f(r)=0$ and $\Delta V(r)=\tilde{V}_{1}-\tilde{V}_{2}$, which is identical to the expression for the conductors discussed above. More general cases are discussed below.

The solution to eq. (4) is uniquely determined by the normalized distance $h/a$ and the average permittivity, $\epsilon_{\rm r}\equiv 2/(\epsilon_{\rm r,1}^{-1}+\epsilon_{\rm r,2}^{-1})$. As shown by Batchelor,¹⁶ eq. (4) can be non-dimensionalized to

$$f(\rho)=\frac{1}{\pi}\int_{0}^{\infty}\!{\rm d}\rho^{\prime}\,\frac{1-f(\rho^{\prime})}{\lambda+\rho^{\prime 2}}\frac{4\rho^{\prime}}{\rho+\rho^{\prime}}K(x)$$

(5)

in which $\rho\equiv r\epsilon_{\rm r}/a$ and $\lambda\equiv h\epsilon_{\rm r}^{2}/a$. The numerically solved $f(\rho)$ for a few representative $\lambda$ values are shown in Fig. 2a. For large gap distance, with $\lambda\gg 1$, $f(\rho)$ is nearly uniform, as expected. For smaller $\lambda$ values, $f(\rho)$ decreases from $f(0)$ with $r$ monotonically. The difference in electric potential at the contact points is proportional to $1-f(0)$. The value of $f(0)$ increases as $h$ decreases, and reaches unity at $h=0$, ensuring that the surface potential is continuous at the contact point. The variation of $f(\rho)$ obtained from the above local analysis is confirmed by directly solving the full potential distribution for dielectric dimers (inset, Fig. 2a).

Similar to eq. (1) for the conductor case, the singular part of surface charges on particle $1$ is given, with $R_{0}$ being the regularizing cutoff of order $a$, by

$$Q_{\rm s,1}=\frac{(\tilde{V}_{1}-\tilde{V}_{2})\epsilon_{\rm out}}{4\pi}\left(1-\epsilon_{{\rm r},1}^{-1}\right)\int_{0}^{R_{0}}\!\!{\rm d}r\frac{2\pi r\left[1-f(r)\right]}{h+r^{2}/a}$$

(6)

The dependence on the dielectric permittivity appears in the prefactor $1-1/\epsilon_{{\rm r},1}$ and in $f(r)$. The singular surface charge $Q_{\rm s,2}$ on the particle 2 is given analogously, but with a negative sign. However, because of the dependence on the dielectric permittivity, $Q_{{\rm s},2}$ and $Q_{{\rm s},1}$ do not add up to zero. The first term in the square bracket gives the same $\ln(a/h)$ singular contributions as eq. (1). The second term represents the correction due to dielectric screening.

The singular capacitance defined by $c_{1}(h)\equiv Q_{{\rm s},1}/(\tilde{V}_{1}-\tilde{V}_{2})$ also has these two contributions. In the non-dimensionalized form, it reads

	$\displaystyle c_{1}(h)$	$\displaystyle=\frac{1-\epsilon_{{\rm r},1}^{-1}}{4}a\epsilon_{\rm out}\left(\ln\frac{a}{h}-P\left(\lambda\right)\right),\quad$		(7)
	$\displaystyle P(\lambda)$	$\displaystyle\equiv\int_{0}^{\infty}\!{\rm d}\rho\,\frac{2\rho}{\lambda+\rho^{2}}f(\rho).$

In the last term, the upper bound is set to infinity because $f(\rho)$ decays rapidly (as $\ln\rho/\rho$, ref.16). The dielectric correction is contained in the term $P(\lambda)$, which depends on permittivity $\epsilon_{\rm r}$ and relative gap distance $h/a$ through the combination $\lambda=h\epsilon_{\rm r}^{2}/a$. The behavior of $P(\lambda)$ is derived from that of $f(\rho)$. For $\lambda\gg 1$, $P(\lambda)$ is vanishingly small, so the capacitance is dominated by $\ln(a/h)$, the characteristic behavior of conductors. On the other hand, when $\lambda\to 0$, the leading contribution to $P(\lambda)$ is $-\ln\lambda$, which cancels exactly the $\ln(a/h)$ dependence, leaving a term $2\ln\epsilon_{\rm r}$ that diverges instead with the average permittivity. Further, we note that $c_{1}(h)$ approaches the conductor limit for $\epsilon_{\rm r}\gg 1$.

The difference in contact charges between dielectric and conductor cases is solely contained in $P(\lambda)$. For intermediate separation, $c_{1}(h)$ shows a singular $\ln(a/h)$ dependence that is similar to the conductor behavior. However, as long as $\epsilon_{\rm r}$ is finite, this logarithmic behavior will eventually be cut off by contribution from $P(\lambda)$ at sufficiently small separation. Therefore, unlike the conductor case, the contact capacitance for dielectrics is finite at $h=0$. Instead, it approaches a constant proportional to $2\ln\epsilon_{\rm r}$. Physically, the logarithmic $h$-dependence originates from the accumulated polarization charges at the interface. It is cut off for dielectrics because the potential difference, which gives rise to the polarization charge, is self-consistently determined by the latter. The weaker polarization effect for dielectrics eventually can not keep up with the needed potential difference for producing the polarization charge. The variation of $c_{1}$ on $h$ and $\epsilon_{\rm r}$ are shown in Fig. 2b. The crossover can be estimated by setting $\lambda=h\epsilon_{\rm r}^{2}/a=1$. The logarithmic divergence, although being cut at small gap separation, still plagues numerical calculations in practice.¹⁸

This type of crossover behavior has been confirmed in a study on dielectric spherical dimers using the image method.¹⁹ The limiting value of the capacitance $2\ln\epsilon_{\rm r}$ as well as the crossover is also consistent with the analytically known result between a conducting sphere and a dielectric plane.²⁰ However, unlike these two work on dimer particles, Batchelor’s result is strictly local, suggesting that the same type of singularity as represented by eq. (7) is applicable to all contact region in cluster of multiple dielectric particles, irrespective of the particle shape. In the following, we show that, the approach for isolating the singularity for contact between conductors can be generalized to treat the dielectric clusters, yielding the exact contact energy and distance-dependence with modest numerical cost.

We consider the cluster consisting of $n$ dielectric spheres. Given the free charge distribution $\rho\lower 3.0pt\hbox{${}_{\rm f}$}({\bf r})$, the potential $\phi$ is governed by Poisson’s equation $\nabla\cdot\epsilon({\bf r})\nabla\phi=-{\rho\lower 3.0pt\hbox{${}_{\rm f}$}({\bf r})}/{\epsilon_{\rm out}}$, where $\epsilon_{\rm out}$ is the medium permittivity. The (relative) material permittivity $\epsilon({\bf r})$ is set to $\epsilon_{\rm r}\equiv\epsilon_{\rm in}/\epsilon_{\rm out}$ inside particles and unity in the medium. Although the general charge distribution poses no added difficulty, for simplicity, we focus on the case when only the free surface charge $\sigma_{\rm f}$ are present. In such case, the system energy can be written as surface integrals of the product of the free surface charge and the surface potential over all surfaces. Furthermore, when the surface charge distribution is uniform, the energy reduces to the sum of the product of the total charge and the average surface potential. Therefore, we denote the set of net charges on particles by ${{\bf Q}}$, and the mean surface potentials ${{\bf V}}$. The net charges and the mean potentials are linearly related, i.e.,

$${\bf Q}={\bf C}{\bf V},$$

(8)

where ${\bf C}$ is an $n{\times}n$ “capacitance array”. In this notation, the total energy can be expressed as $\mathcal{E}=\frac{1}{2}{{\bf Q}}\cdot{\bf C}^{-1}\cdot{\bf Q}$. For convenience, the total energy $\mathcal{E}$ presented below is always normalized by the self-energy of a single isolated sphere with the radius $a$ and a net charge $q$, $q^{2}/(8\pi\epsilon_{\rm out}a)$. We note that this formulation applies to arbitrary particle shapes and cluster configurations. If the surface charge distribution is non-uniform or any external excitation exists, a multipole expansion of surface charge in terms of dipole, quadrupole etc. is needed. The constitutive relation eq. (8) and the quadratic expansion to energy can be generalized to include the contribution from these higher multipoles and external electric fields. For the purpose of demonstrating how the contact singularity is isolated, we focus on the case of the uniform surface charge distribution.

Our method for dielectrics is an extension to our earlier work on conductors.⁵ Because the conductors are equipotential, the mean potential is also the surface potential at every point on the surface. Since the capacitance $\bf C$ in eq. (8) contains two types of contributions, one type dominated by close contacts between neighboring particles, and the other type from the remaining interactions from all particles, we decompose the capacitance $\bf C$ as follows

$${\bf Q}=\left(c_{\rm s}{\bf L}+{\bf H}\right){\bf V}.$$

(9)

Here, $c_{\rm s}(h)=\frac{a\epsilon_{\rm out}}{4}\ln(a/h)$ is the singular capacitance for conductors derived above (assuming that all gap distances are $h$). Since $c_{\rm s}$ becomes significant only for small gaps, the entries in the array $\bf L$ are nonvanishing only for closely neighboring particles. Specifically, if particle $i$ and $j$ form a close contact, we set the entries $L_{ij}=L_{ji}=-1$. The diagonal entry $L_{ii}$ equals the number of close neighbors of the particle $i$. All other entries of $\bf L$ are set to zero. The singular capacitance $\bf L$ as defined here is the same as the adjacency matrix representing the connectivity of clusters (see ref.5 for explicit examples). In practice, for a given cluster configuration, we first numerically solve the full capacitance array $\bf C$ at finite but small $h$ values, then subtract from $\bf C$ the singular term $c_{\rm s}(h){\bf L}$, to get the regular capacitance ${\bf H}$. The $h$-dependence of this regular capacitance is then fitted to a straight line when $h$ is small, allowing us to obtain the full $h$-dependence for the capacitance array $\bf C$ and, consequently, the full $h$-dependence of energy, down to $h=0$.

Generalizing eq. (9) to dielectrics requires two nontrivial modifications. First, the decomposition in eq. (9) is valid because the potential of the conductors can be used to evaluate the contact potential difference. But dielectrics are not equipotential, and the surface potential at the contact points $\tilde{V}_{i}$ generally differ from the average potential $V_{i}$ by a numerical factor that depends on the dielectric permittivity and the cluster configuration. Its variation with gap distance $h$ is weak, and approaches a constant as $h\to 0$. So we generalize eq. (9) to

$${\bf Q}=\left(\hbox{$\bf\tilde{\phantom{L}}$\kern-5.50003pt\hbox{$\bf L$}}+{\bf H}\right){\bf V},\quad\tilde{L}_{ij}\equiv\frac{1-\epsilon_{{\rm r},j}^{-1}}{V_{j}}\begin{cases}\tilde{V}_{j}^{(i)}c_{ij},&i\neq j;\\[5.0pt] {\sum_{k=1}^{N}}\tilde{V}_{j}^{(k)}c_{jk},&i=j.\end{cases}$$

(10)

Above, the superscript ‘$(i)$’ in $\tilde{V}_{j}^{(i)}$ indicates that the contact potential is evaluated on the particle $j$ at the contact formed with the particle $i$. The singular capacitance $c_{ij}$ is given from eq. (7) by $c_{ij}=\frac{\epsilon_{\rm out}a_{ij}}{4}\left(\ln\frac{a_{ij}}{h_{ij}}-P(\lambda_{ij})\right)$, which depends on the mean radius of curvature $a$, gap distance $h$, and $\lambda$ value evaluated for the particle pair $i$ and $j$. From the definition, it is clear that $\bf\tilde{\phantom{L}}$$\bf L$ is generally not symmetric, which however reduces to the (symmetric) form identical to eq. (9) in the conductor limit, since $\tilde{V}_{j}^{(i)}=V_{j}$ for all contacts on the particle $j$.

Second, it turns out that the contact potentials exhibit strong dependence on the second order and more distant neighbors, because the dielectric screening is comparatively weak (see Fig. 4 below). Therefore, to ensure rapid convergence of the correct contact potentials, our construction of the adjacency array $\bf\tilde{\phantom{L}}$$\bf L$ contains all pairs of particles. Fortunately, as we shall see below, this construction requires no extra computation other than evaluating the surface potentials at additional contact points.

In order to identify the regular part ${\bf H}$ using eq. (10), we numerically compute the full capacitance matrix ${\bf C}(h)$ and all the contact potentials $\tilde{V}_{j}^{(i)}$, for a range of small but nonzero $h$ values. We used the boundary element method (BEM) implemented in the package COPSS ¹⁸ to solve Poisson’s equation. The solution is expressed in terms of the induced bound (polarization) charges $\sigma_{\rm b}$, which satisfies the following boundary integral equation,

$$\frac{1+\epsilon_{\rm r}}{2}\sigma_{\rm b}+(1-\epsilon_{\rm r})\epsilon_{\rm out}\left({\bf E}_{\rm b}+{\bf E}_{\rm f}\right)\cdot\hat{\bf n}=\frac{1-\epsilon_{\rm r}}{2}\sigma_{\rm f}.$$

(11)

Above, the electric field strength at the surface due to induced bound charge $\sigma_{\rm b}$ and free charges $\sigma_{\rm f}$ are expressed as surface integrals ($\alpha={\rm b,f}$),

$${\bf E}_{\alpha}=\frac{1}{4\pi\epsilon_{\rm out}}\int\!{\rm d}S^{\prime}\,\frac{{\bf r}-{\bf r}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|^{3}}\sigma_{\alpha}({\bf r}^{\prime}).$$

(12)

By our convention, the surface normal $\hat{\bf n}$ points outward. For clusters of multiple particles, it is understood that different particles may have different values of permittivity $\epsilon_{\rm r}$. Equation (11) is linear in $\sigma_{\rm b}$ and $\sigma_{\rm f}$. After discretization, the bound charge is obtained by inverting the coefficient array. The apparent divergent diagonal entries while discretizing the surface integral eq. (12) can be re-parameterized to reproduce the correct self-energy. Further numerical details can be found in ref.8. In all calculations by BEM presented in this work (the dimer example used the exact series expansion.¹⁰), we meshed each spherical surface into $17342$ triangular patches such that mesh size is about $0.03\,a$.

In general, for a cluster of $n$ particles, $n$ separate numerical calculations with $n$ independent charge vectors $\bf{Q}$ are needed to determine the full capacitance matrix ${\bf C}$. The computational cost can be reduced by imposing the symmetry of the cluster configuration. Subtracting from $\bf C$ the singular contribution $\bf\tilde{\phantom{L}}$$\bf L$ is expected to give the $h$-dependence of the regular part $\bf H$. However, unlike the conductor case, where the coefficient to the logarithmic term is exactly known, the singular capacitance $\bf\tilde{\phantom{L}}$$\bf L$ for dielectrics depends on the calculated contact potentials $\tilde{V}_{j}^{(i)}$, which is susceptible to the numerical precision and what is meant by ‘contact point’ on a surface mesh. In practice, we vary the gap distance $h$ and compute a few trial values of contact potentials, then select the one that ensures the difference $H_{ij}(h)=C_{ij}-\tilde{L}_{ij}$ converges linearly with $h$ as $h\to 0$ for all the entries. In the following, three examples are presented to demonstrate the nature of the contact singularity, and to show that the regularization scheme allows us to obtain the energy of cluster of dielectric particles with small $h$.

The first example is a dimer of identical dielectric spheres. It is the simplest example for demonstrating the effects of interfacial polarization.^{6, 9, 10, 21} Except a study on the interaction between a conducting sphere and a dielectric plane,²⁰ very few past work tried to evaluate the energy at small separation, presumably due to the difficulty of resolving the aforementioned contact singularity. Therefore, we analyzed the polarization effect for dielectric dimers in close contact, verified the singular behavior in eq. (7), and showed that eq. (10) yields the full $h$-dependence of energy. We studied both symmetric (${\bf Q}=(q,q)$) and asymmetric (${\bf Q}=(q,-q)$) cases, for a range of permittivity values ($1\leq\epsilon_{\rm r}\leq 10^{6}$). The energy at $h=0$ is presented as a function of relative permittivity, which agrees with the exact result (see SI) found using the tangent-sphere coordinate.

For the asymmetric case, interfacial polarization enhances the inter-particle attraction. To illustrate this effect, we subtract the energy of two isolated spheres from that of a dimer, then plot it against the gap distance in Fig. 3, for three representative values of $\epsilon_{\rm r}$: $10$, $10^{2}$ and $10^{4}$. In absence of the polarization effect, the energy curves would exhibit no dependence on $\epsilon_{\rm r}$, and appear flat over this narrow range of $h/a$ values. As $\epsilon_{\rm r}$ is increased from $1$, we confirm the expected, stronger distance dependence, as $h\to 0$. In particular, for the case of $\epsilon_{\rm r}=10^{4}$, an apparent logarithmic dependence on $h$ is seen in Fig. 3a, which is precisely the singularity revealed in eq. (7) and is stronger for higher permittivity values. On the other hand, when the permittivity $\epsilon_{\rm r}$ is decreased, this logarithmic dependence is cut off at an increasingly larger separation and becomes almost invisible when $\epsilon_{\rm r}=10$.

The normalized energy difference, as shown in Fig. 3a, evaluated at $h=0$, is show in Fig. 3b for different $\epsilon_{\rm r}$ values. This contact energy represents the energy gain for bringing two particles from infinity to contact. In terms of the normalizing energy unit, $q^{2}/(8\pi\epsilon_{\rm out}a)$, its value is $-1$ without the polarization effect, and becomes more negative as $\epsilon_{\rm r}$ increases, eventually reaching $-2$ at $\epsilon_{\rm r}=\infty$.⁵ The contribution from the polarization effect for $\epsilon_{\rm r}\gtrsim 100$ is comparable to the coulombic attraction alone. Moreover, we notice a slow convergence of contact energy for permittivities $\epsilon_{\rm r}\gtrsim 10^{3}$, owing to the weak $\ln\epsilon_{\rm r}^{2}$ dependence in the contact capacitance eq. (7) at $h=0$. Finally, we affirm our results obtained from eq. (10), by comparing the contact energies to the exact analytical results (SI).

For the symmetric case, the interfacial polarization weakens the inter-particle repulsion, which is seen from the $h$-dependence of dimer energy. However, by symmetry, the potential difference in the gap region vanishes and thus no logarithmic behavior is observed in Fig. 3c. The energy scales linearly with $h$ and a straightforward extrapolation gives the contact energy for all permittivity values plotted in Fig. 3d. A much weaker polarization effect is found. Even for the conductor case, only about $10\%$ contact energy can be attributed to polarization effects. The contact energies converge for $\epsilon_{\rm r}\gtrsim 100$, and reach $2(1/\ln 2-1)\simeq 0.88$, in agreement with Maxwell’s result²² for the conducting dimers. As for the asymmetric case, the full $\epsilon_{\rm r}$-dependence matches the analytical calculation (SI). To conclude this example, we note that the symmetric case is a peculiar example where the contact singularity is strictly absent. For all other charge ratios, using the decomposition in eq. (10) to isolate the strong $h$-dependence is necessary.

The second example concerns the energy of $8$ identical spheres placed at the vertices of a cube, which illustrates the importance of the second order contact for dielectric clusters. As discussed above, the singular capacitance $\bf\tilde{\phantom{L}}$$\bf L$ in eq. (10) contains not only contributions from the nearest neighbors, but also those from the second-order and higher order neighbors. Therefore, a sphere at a vertex of a cube forms a secondary contact with its $3$ plane diagonal vertices, and a third order contact with the body diagonal vertex. Even though the higher order contributions are minor, because at such large separation, these contacts cease to be singular, keeping the second order contribution is essential for obtaining the correct linear scaling of the regular capacitance with gap distance $h$ shown in Fig. 4a.

As in the dimer case, these regularized capacitance allow us to evaluate the energy at arbitrarily small gap distance. Figure 4b shows the normalized energies when only one particle is charged, for which all the energetic contributions come from the interfacial polarization. Comparing the results from BEM and our regularization schemes containing varying order of contact contributions, it is clear that keeping the higher order singular contribution is necessary for correctly evaluating the energy for $h/a\lesssim 0.05$. In contrast, no such terms are needed for the conducting case.

The third example is our main result on the energy of clusters, from which the cohesive energy is obtained by subtracting the self-energy of particles at infinite separation. We considered two limiting configurations: the most extended one with all spheres arranged along a straight line (string), and the most compact one with all spheres posited at the vertices of the platonic solids (polyhedron). In our earlier work on the conducting spheres, the charges are allowed to flow freely between contacting spheres. The energy of polyhedron packing is found to be lower than that of the string packing at a finite separation. However, as the gap distance decreases, the polarization effects due to contact singularity become increasingly relevant, which ultimately makes the string packing to be energetically more favorable than the polyhedron packing. For the dielectric cases, we show that the weakened contact singularity causes another crossover between the relative stability of string-like and polyhedral packings.

As the previous example, to highlight the polarization effects, a unit charge $q$ is placed on one sphere. For symmetric polyhedron packing, this charged sphere is arbitrarily chosen. For the string packing, we considered two extreme placings: at one end (string-1) and in the middle (string-2). For each particle configuration and charge assignment, we considered two scenarios: charge transfer between contacting particles is permitted (Fig. 5a) and prohibited (Fig. 5b).

When inter-particle charge transfer is permitted, while maintaining uniform charge distribution on each individual particle, the energy $\mathcal{E}$ can be calculated using ${\bf Q}\cdot{\bf C}^{-1}\cdot{\bf Q}/2$. Minimizing this energy subject to the constraint of constant total charge $q$, the optimal charge assignment is found to be ${\bf Q}_{e}=q\,{\bf C}\cdot{\bf u}/C_{\rm e}$, where ${\bf u}\equiv(1,1,\cdots,1)$ and $C_{\rm e}\equiv\sum_{i,j=1}^{n}C_{ij}$. This implies that the optimal charge ${\bf Q}_{e}$ produces identical average surface potentials for all the particles, which generalizes the ‘equipotential’ concept for conductors.⁵ Correspondingly, the minimized energy is given by $q^{2}/(2C_{\rm e})$, where the mean capacitance $C_{\rm e}$ weakly depends on the dielectric permittivity. In this scenario, there is no need to differentiate string-$1$ and string-$2$ configurations. The energies of string and polyhedron configurations at $h=0$ are plotted in Fig. 5a, for  $\epsilon_{\rm r}=10$ and $500$. For all configurations, the energy decreases with the cluster size $n$ monotonically since the redistributed surface charges are further apart. The dependence on permittivity is surprisingly weak, and the size-dependence is nearly indistinguishable from that of the conductor case.⁵ Two curves are obtained by treating the cluster as a single conductor respectively, whose capacitance scales with the length scale of the cluster, i.e., $n^{1/3}$ for the polyhedron packing and $n/\ln(n)$ for the string packing.⁵ The extended string packing has a lower cohesive energy because its effective capacitance is higher than that of the compact polyhedron packing.

More interesting behavior is found (Fig. 5b) when charge transfer is prohibited. The dependence on permittivity is evident for all three cases: string-$1$, string-$2$, and polyhedron. To better visualize energies of string-$2$ and polyhedron packing, energies of string-$1$ is presented in Fig. S4) The energy is lower for the cluster with higher permittivity, because the polarization effect is stronger. The energy of string-$2$ is lower than that of string-$1$ for identical $n$, because the charged particles in the middle of the string can polarize the particles in two half strings. Further, the energies for both string-$1$ and string-$2$ configurations saturate as the cluster size grows beyond $n=12$, because the polarization effect is short-ranged. To assess the relative stability of compact or extended configurations, we then only need to focus on the string-$2$ and the polyhedron cases.

Our results indicate that the relative energetic stability depends on the permittivity. For $\epsilon_{\rm r}=10$, the energy of the polyhedron packing is lower than the string-2 packing, which is opposite to the characteristic result for conductors shown in Fig. 5a. Therefore, we expect a crossover from stable compact packing to stable open packing at intermediate permittivity values. This is indeed the case shown for $\epsilon_{\rm r}=500$, whereby the energies of compact and open configurations closely trace each other, and the relative energetic stability changes at $n=4$ and $n=8$ respectively. As $\epsilon_{\rm r}$ is further increased, the energy curves of corresponding configuration will converge to those in Fig. 5a, reversing the relative stability of open and compact packings.

In summary, we generalized our previous work on conducting particles,⁵ and developed a scheme to resolve the singular contact charges between touching dielectric spheres, which regularizes the full capacitance array by isolating the singular contributions, i.e., eq. (10). Using this scheme, we obtain the cohesive energies for dielectric clusters at zero separation containing up to $n=20$ particles, which is difficult to resolve with brutal force numerical calculations. Our results show that the shape of stable clusters formed from dielectric particles depends on the permittivity ratio $\epsilon_{\rm r}$: open clusters is more stable for large $\epsilon_{\rm r}$, and compact clusters is more stable for small $\epsilon_{\rm r}$. Our scheme is applicable to systems with arbitrary packing geometry, charge distribution, and containing asymmetric interfaces that have different permittivities and radii of curvature on the contacting particles. At the heart of our scheme is the contact singularity, which resembles those encountered in the study of thermal conduction ¹⁶ and momentum transport across a narrow gap.²⁴