Debye-Hückel Free Energy of an Electric Double Layer with Discrete Charges Located at a Dielectric Interface

Abstract

Poisson–Boltzmann theory provides an established framework to calculate properties and free energies of an electric double layer, especially for simple geometries and interfaces that carry continuous charge densities. At sufficiently small length scales, however, the discreteness of the surface charges cannot be neglected. We consider a planar dielectric interface that separates a salt-containing aqueous phase from a medium of low dielectric constant and carries discrete surface charges of fixed density. Within the linear Debye-Hückel limit of Poisson–Boltzmann theory, we calculate the surface potential inside a Wigner–Seitz cell that is produced by all surface charges outside the cell using a Fourier-Bessel series and a Hankel transformation. From the surface potential, we obtain the Debye-Hückel free energy of the electric double layer, which we compare with the corresponding expression in the continuum limit. Differences arise for sufficiently small charge densities, where we show that the dominating interaction is dipolar, arising from the dipoles formed by the surface charges and associated counterions. This interaction propagates through the medium of a low dielectric constant and alters the continuum power of two dependence of the free energy on the surface charge density to a power of 2.5 law.

1. Introduction

Charged interfaces in aqueous solution give rise to the formation of an electric double layer—a diffuse cloud of co- and counter-ions that screen the interfacial charges [1,2]. Electric double layers are abundant in every living cell, as well as in many technological applications. Among the most prominent examples are biomembranes, which consist of a mixture of lipids and associated proteins that self-assemble into an extended planar sheet. Some of the lipids are usually anionic, with their headgroup charges being located at the dielectric interface between a region of low dielectric constant inside the membrane and the aqueous phase, which has a large dielectric constant and contains salt ions [3].

The properties of the electric double layer can be described at different levels of abstraction, varying from detailed atomistic and coarse-grained computer simulations [4] to the most simple mean-field models. A major feature of the latter is the neglect of ion-ion correlations [5]. Mean-field models become especially useful [6,7] and amenable to extensions [8,9] when based on simple geometries and when treating the dielectric interface as uniformly charged. One of the most important extensions of applying mean-field electrostatics to the electric double layer involves inhomogeneously charged surfaces, ranging from continuous charge distributions to discrete point charges. Surfaces with patchy or non-uniform charge distributions [10,11,12,13], domain formation and phase separations [14,15], and domain boundaries [16,17,18] all have been modeled using mean-field electrostatics while maintaining the assumption of charge continuity, which is justified on sufficiently large length scales. For weakly charged surfaces or on sufficiently small length scales, the discreteness of the electric charge will become noticeable. This has inspired a number of authors to compute the electrostatic potential produced by an ordered array of surface charges [19,20,21] or even a single surface charge [22,23,24,25] and others to use virial expansions [26,27] and Monte Carlo simulations [28,29,30,31,32,33,34,35,36,37] to account for charge discreteness. Some studies have focused on the pressure that acts across electrolytes due to the presence of discrete charges [28,35,38], yet—perhaps somewhat surprisingly—there has not been an attempt so far to compute the free energy of a single planar dielectric interface with discrete charges and compare this with the corresponding free energy derived for a continuous charge distribution. Filling this gap on the level of the linearized Debye-Hückel model is the goal of the present work.

We study the properties of an electric double layer on the basis of the linearized Debye-Hückel limit, which becomes valid for a sufficiently small average surface charge density. Discrete charges are attached with a certain density to a planar dielectric interface that separates a salt-containing aqueous solution from a medium with a low dielectric constant. We introduce a Wigner–Seitz cell model [39] to account for the interactions between the interfacial charges, which do not need to form an ordered array. Using a Fourier-Bessel series and a Hankel transformation, we calculate the electrostatic surface potential at the position of an individual charge produced by all other charges on the interface and, from that, the Debye-Hückel free energy of the electric double layer. The discrete nature of the interfacial charges becomes important for sufficiently small densities. In this case, the interactions between the interfacial charges are dominated by dipolar contributions mediated through the salt-free medium of low dielectric constant.

2. Theory and Discussion

We consider a flat interface of sufficiently large lateral area A that separates two distinct dielectric media. One medium has a dielectric constant ${\epsilon }_l$, and the other medium has a dielectric constant ${\epsilon }_w$ and contains monovalent anions and cations of a salt such as NaCl, both of bulk concentration $n_0$. Because a lipid layer is representative of the systems we intend to model, we chose the subscripts “l” and “w” for the “lipid” hydrocarbon core and “water”, respectively, where $2\lesssim {\epsilon }_l\lesssim 4$ and ${\epsilon }_w\approx 80$. However, our theoretical model is general and applies to any choices of ${\epsilon }_l$ and ${\epsilon }_w$. The presence of salt ions entails the formation of an electric double layer in the medium of dielectric constant ${\epsilon }_w$. The characteristic screening length of the electric double layer is given by the Debye length $l_D={[\left({\epsilon }_w{\epsilon }_0{k}_{B}T\right)/\left(2e^2{n}_0\right)]}^{1/2}$, where ${\epsilon }_0$ is the vacuum permittivity, $k_B$ Boltzmann’s constant, T the absolute temperature, and e the elementary charge.

The dielectric interface also carries N uniformly distributed elementary charges e, either positive or negative ones, implying an average surface charge density ${\sigma }_0=\pm eN/A$. The charges exhibit some short-range order (as is the case in a fluid); long-range order is not necessary for our model to apply. For a lipid layer, $|{\sigma }_0|\lesssim e/{nm}^2$, with the exact value depending on the fraction of charged lipids. When ${\sigma }_0$ is sufficiently small, linearized Poisson–Boltzmann theory—also known as the Debye-Hückel limit—offers a simple and thus commonly used model to describe properties of the electric double layer such as the electrostatic surface potential ${\mathsf{\Phi }}_0={\sigma }_0{l}_D/\left({\epsilon }_0{\epsilon }_w\right)$ or the free energy per unit area:

$$\frac{F_0}{A}=\frac{{\sigma }_0^2}{2}\frac{l_D}{{\epsilon }_w{\epsilon }_0}$$

(1)

of the charged interface. The validity of Equation (1) rests on a number of assumptions, including the absence of ion-ion correlations, a uniform dielectric background (of dielectric constant ${\epsilon }_w$), and a sufficiently small and yet uniform and continuous charge density ${\sigma }_0$ at the dielectric interface. The latter assumption must break down at some point due to the discrete nature of the charges at the interface. That is, when the separation between individual charges is much larger than the Debye screening length, these charges do not interact anymore, suggesting that Equation (1) overestimates the free energy $F/A$ that would account properly for charge discreteness. In the following, we calculate the free energy $F/A$ for a set of uniformly distributed discrete charges at the dielectric interface and compare it with $F_0/A$ in Equation (1), which ignores the discrete nature of the charges.

2.1. Ensemble of Charged Disks

Our method to account for the discreteness of the charges on the interface invokes the consideration of Wigner–Seitz cells [40]. A Wigner–Seitz cell serves as a unit cell of cylindrical symmetry with a circular cross-section of radius R and corresponding lateral area $A/N=\pi R^2$. We identify the z-axis of a cylindrical coordinate system $(r,\varphi ,z)$ with the symmetry axis of the Wigner–Seitz cell and locate the plane $z\equiv 0$ at the dielectric interface. The radial coordinate r varies in the region $0\le r\le R$, from the symmetry axis to the boundary of the unit cell. Note that the cylindrical symmetry of the unit cell renders all system properties invariant with respect to the azimuthal angle $\varphi$.

If the charge on the dielectric interface was continuous and distributed uniformly, the surface charge density $\sigma \left(r\right)={\sigma }_0$ along the radial direction would be strictly constant. In contrast, we shall consider a non-uniform, step-like surface charge density:

$$\sigma \left(r\right)=\{\begin{array}{ccc}\overline{\sigma }\hfill & \mathrm{if}& r\le r_0\hfill \\ 0\hfill & \mathrm{if}& r_0<r\le R,\hfill \end{array}$$

(2)

with $\overline{\sigma }={\sigma }_0{R}^2/r_0^2$. That is, we condense all charge, ${\sigma }_0\pi R^2$, contained in the unit cell into a concentric circular region of radius $r_0$, leaving the outer region $r_0<r\le R$ of the dielectric interface inside the unit cell uncharged. We are free to choose both $r_0$ and the total amount of charge in the unit cell ${\sigma }_0\pi R^2$. Yet, in order to study the influence of charge discreteness, we will eventually identify ${\sigma }_0\pi R^2=\overline{\sigma }\pi r_0^2$ with a single elementary charge e and consider the limit $r_0\to 0$ of a point charge. Our goal is to calculate the corresponding Debye-Hückel free energy $F\left(r_0\right)$ in that limit. The left side of Figure 1 shows an illustration of the dielectric interface (red plane) with seven circular disks of radius $r_0$ (shown in blue color) that carry a charge density $\overline{\sigma }$ each, one with its unit cell (the Wigner–Seitz cell, marked by the blue transparent cylinder) explicitly shown.

The right side of Figure 1 presents the cylindrical coordinate system associated with the unit cell; the outer boundary, at $r=R$, is marked by a dashed line, and the interface is located at the plane $z\equiv 0$.

We denote the electrostatic potential for $z\ge 0$ by ${\mathsf{\Phi }}_1(r,z)$. Due to the presence of salt, ${\mathsf{\Phi }}_1$ fulfills the Debye-Hückel equation $l_D^2{\nabla }^2{\mathsf{\Phi }}_1={\mathsf{\Phi }}_1$, where ${\nabla }^2$ denotes the Laplacian operator. Similarly, we denote the electrostatic potential for $z\le 0$ by ${\mathsf{\Phi }}_2(r,z)$. The potential ${\mathsf{\Phi }}_2$ fulfills the Laplace equation ${\nabla }^2{\mathsf{\Phi }}_2=0$. The Debye-Hückel and Laplace equations for a cylindrically symmetric unit cell and the structure of their solution can be written as:

$$\begin{array}{ccc}\hfill \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial {\mathsf{\Phi }}_1}{\partial r}\right)+\frac{{\partial }^2{\mathsf{\Phi }}_1}{\partial z^2}=\frac{{\mathsf{\Phi }}_1}{l_D^2}& \to & {\mathsf{\Phi }}_1(r,z)=\sum _{n=0}^{\infty }{a}_n{J}_0\left(y_n\frac{r}{R}\right)e^{-z\sqrt{{\left(\frac{y_n}{R}\right)}^2+\frac{1}{l_D^2}}}\hfill \\ \hfill \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial {\mathsf{\Phi }}_2}{\partial r}\right)+\frac{{\partial }^2{\mathsf{\Phi }}_2}{\partial z^2}=0& \to & {\mathsf{\Phi }}_2(r,z)=\sum _{n=0}^{\infty }{a}_n{J}_0\left(y_n\frac{r}{R}\right)e^{z\frac{y_n}{R}},\hfill \end{array}$$

(3)

where $J_0$ denotes the Bessel function of the first kind and zeroth order. The sets of the yet undetermined constants $a_n$ in the expressions for ${\mathsf{\Phi }}_1(r,z)$ and ${\mathsf{\Phi }}_2(r,z)$ are identical to ensure the electrostatic potential remains continuous when passing through the dielectric interface at $z=0$. In addition, the solutions do not diverge for $\left|z\right|\to \infty$. Symmetry requires the derivatives ${(\partial {\mathsf{\Phi }}_1/\partial r)}_{r=R}$ and ${(\partial {\mathsf{\Phi }}_2/\partial r)}_{r=R}$ at the boundary of the unit cell to vanish. Because the Bessel function of the first kind and first order, $J_1\left(y\right)=-J_0^{\prime }\left(y\right)$, is the negative derivative of the Bessel function of the first kind and zeroth order, we demand $J_1\left(y_n\right)$ to vanish for all $y_n$ with integers $n\ge 0$. Hence, $y_n$ in Equation (3) satisfies $J_1\left(y_n\right)=0$ and is thus the n’th zero of the Bessel function of the first kind and first order. Figure 2 shows $J_0\left(y\right)$ and $J_1\left(y\right)$, with the set $y_n$ indicated on the upper axis for $0\le n\le 9$.

At the dielectric interface, the electrostatic potential must fulfill the boundary condition:

$${\epsilon }_w{\left(\frac{\partial {\mathsf{\Phi }}_1}{\partial z}\right)}_{z=0}-{\epsilon }_l{\left(\frac{\partial {\mathsf{\Phi }}_2}{\partial z}\right)}_{z=0}=-\frac{\sigma \left(r\right)}{{\epsilon }_0}.$$

(4)

The functions $J_0(y_{n}r/R)$ can be used to represent a given surface charge density $\sigma \left(r\right)$ through a Fourier-Bessel series $\sigma \left(r\right)={\sum }_{n=0}^{\infty }{\sigma }_n{J}_0(y_{n}r/R)$ in the region $0\le r\le R$, with a set of constants ${\sigma }_n$. Inserting that series and the expressions for ${\mathsf{\Phi }}_1(r,z)$ and ${\mathsf{\Phi }}_2(r,z)$ from Equations (3) into the boundary condition specified by Equation (4), yields the relation:

$$a_n=\frac{1}{{\epsilon }_0}\frac{{\sigma }_n}{{\epsilon }_w\sqrt{{\left(\frac{y_n}{R}\right)}^2+\frac{1}{l_D^2}}+{\epsilon }_l\frac{y_n}{R}}$$

(5)

between the coefficients $a_n$ and ${\sigma }_n$ for all integers $n\ge 0$. Note that, as required by symmetry, the derivative ${\sigma }^{\prime }\left(r\right)=-{\sum }_{n=1}^{\infty }{\sigma }_n{J}_1(y_{n}r/R)y_n/R=0$ indeed vanishes at the boundary of the unit cell. Because of ${\int }_0^{R}drr\sigma \left(r\right)={\sum }_{n=0}^{\infty }{\sigma }_n{J}_1\left(y_n\right)R^2/y_n=R^2{\sigma }_0/2$, we identify the first coefficient:

$${\sigma }_0=\frac{\underset{0}{\overset{R}{\int }}drr\sigma \left(r\right)}{\underset{0}{\overset{R}{\int }}drr}$$

(6)

of the Fourier-Bessel series with the average surface charge density, in agreement with how we have already used it in Equation (1). To find the remaining coefficients ${\sigma }_n$ for $n\ge 1$, we employ the orthogonality ${\int }_0^{R}drr{J}_1(y_{n}r/R)J_1(y_{m}r/R)={\delta }_{mn}{\left[R{J}_0\left(y_n\right)\right]}^2/2$ of the functions $J_1(y_{n}r/R)$ in the region $0\le r\le R$, where ${\delta }_{mn}$ denotes the Kronecker delta. Applying orthogonality to evaluate the integral ${\int }_0^{R}drr{\sigma }^{\prime }\left(r\right)J_1(y_{n}r/R)$ leads to:

$${\sigma }_n=-\frac{\underset{0}{\overset{R}{\int }}drr{\sigma }^{\prime }\left(r\right)J_1\left(y_n\frac{r}{R}\right)}{y_n\frac{R}{2}{\left[J_0\left(y_n\right)\right]}^2}.$$

(7)

Using the results for $a_n$ in Equation (5) and for ${\sigma }_n$ in Equations (6) and (7) gives rise to the explicit expression:

$${\mathsf{\Phi }}_s\left(r\right)={\mathsf{\Phi }}_0-\frac{l_D}{{\epsilon }_0{\epsilon }_w}\sum _{n=1}^{\infty }\frac{J_0\left(y_n\frac{r}{R}\right)}{\sqrt{{\left(y_n\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}{y}_n}\frac{\underset{0}{\overset{R}{\int }}drr{\sigma }^{\prime }\left(r\right)J_1\left(y_n\frac{r}{R}\right)}{y_n\frac{R}{2}{\left[J_0\left(y_n\right)\right]}^2}$$

(8)

for the surface potential ${\mathsf{\Phi }}_s\left(r\right)={\mathsf{\Phi }}_1(r,z=0)={\mathsf{\Phi }}_2(r,z=0)$, valid for any choice of the surface charge density $\sigma \left(r\right)$. Recall that ${\mathsf{\Phi }}_0={\sigma }_0{l}_D/\left({\epsilon }_0{\epsilon }_w\right)$ is the surface potential for a strictly uniform and continuous surface charge density ${\sigma }_0$.

Our next step is to apply the general formalism of calculating the surface potential ${\mathsf{\Phi }}_s\left(r\right)$ in Equation (8) to our specific choice of $\sigma \left(r\right)$ in Equation (2). It is convenient to re-express $\sigma \left(r\right)=\overline{\sigma }u(r_0-r)$ in terms of the Heaviside step function $u\left(r\right)=0$ for $r<0$, $u\left(r\right)=1$ for $r>0$, and $u\left(r\right)=1/2$ for $r=0$. This implies ${\sigma }^{\prime }\left(r\right)=-\overline{\sigma }\delta (r-r_0)$, where $\delta \left(r\right)$ denotes the Dirac delta function. Hence, ${\int }_0^{R}drr{\sigma }^{\prime }\left(r\right)J_1(y_{n}r/R)=-\overline{\sigma }{r}_0{J}_1(y_n{r}_0/R)$, and thus:

$$\begin{array}{ccc}\hfill \frac{\sigma \left(r\right)}{{\sigma }_0}& =& 1+\sum _{n=1}^{\infty }\frac{J_0\left(y_n\frac{r}{R}\right)2\frac{R}{r_0}{J}_1\left(y_n\frac{r_0}{R}\right)}{y_n{\left[J_0\left(y_n\right)\right]}^2},\hfill \\ \hfill \frac{{\mathsf{\Phi }}_s\left(r\right)}{{\mathsf{\Phi }}_0}& =& 1+\sum _{n=1}^{\infty }\frac{J_0\left(y_n\frac{r}{R}\right)}{\sqrt{{\left(y_n\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}{y}_n}\frac{2\frac{R}{r_0}{J}_1\left(y_n\frac{r_0}{R}\right)}{y_n{\left[J_0\left(y_n\right)\right]}^2}.\hfill \end{array}$$

(9)

Equations (9) represent the Fourier-Bessel series of the surface charge density according to Equation (2) as well as the corresponding electrostatic surface potential ${\mathsf{\Phi }}_s\left(r\right)$ at the dielectric interface. Recall that our choice $\overline{\sigma }={\sigma }_0{R}^2/r_0^2$ ensures the total charge in the unit cell is $\pi R^2{\sigma }_0$. Both $\sigma \left(r\right)/{\sigma }_0$ and ${\mathsf{\Phi }}_s\left(r\right)/{\mathsf{\Phi }}_0$ are functions of $r/R$ and depend on the relative size of the charged circular region as compared to the cell size, $r_0/R$. The scaled potential depends on two additional parameters, $l_D/R$ and ${\epsilon }_l/{\epsilon }_w$. Figure 3 shows $\sigma \left(r\right)/{\sigma }_0$ (left diagram) and ${\mathsf{\Phi }}_s\left(r\right)/{\mathsf{\Phi }}_0$ (right diagram) for different choices of $r_0/R$, all with $l_D=R$ and ${\epsilon }_l={\epsilon }_w$.

In the limit $r_0\to 0$, the charge distribution $\sigma \left(r\right)={\sigma }_0\pi R^2\delta \left(r\right)$ approaches that of a single point charge, and the corresponding surface potential becomes:

$$\frac{{\mathsf{\Phi }}_s\left(r\right)}{{\mathsf{\Phi }}_0}=1+\sum _{n=1}^{\infty }\frac{J_0\left(y_n\frac{r}{R}\right)}{\sqrt{{\left(y_n\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}{y}_n}\frac{1}{{\left[J_0\left(y_n\right)\right]}^2}.$$

(10)

Note that ${\mathsf{\Phi }}_s\left(r\right)$ in Equation (10) diverges at the location of the point charge, $r=0$.

2.2. Single Isolated Charged Disk

In order to calculate the interaction energy between an ensemble of charged disks (and eventually an ensemble of point charges), we will also consider a single isolated charged disk as the reference state. We denote the electrostatic surface potential produced by an isolated charged disk of radius $r_0$ by ${\mathsf{\Phi }}_s^{self}\left(r\right)$. The calculation of ${\mathsf{\Phi }}_s^{self}\left(r\right)$ proceeds analogously to that of ${\mathsf{\Phi }}_s\left(r\right)$, yet for an infinitely large unit cell size, $R\to \infty$, while preserving $\overline{\sigma }={\sigma }_0{R}^2/r_0^2$ as specified in Equation (2). In this case, we shall make use of a Hankel transform (also known as Bessel–Fourier transform) instead of a Bessel-Fourier series. For the surface charge density, we employ the identity:

$$\sigma \left(r\right)=\underset{0}{\overset{\infty }{\int }}dkk{J}_0\left(kr\right)\underset{0}{\overset{\infty }{\int }}d\overline{r}\overline{r}\sigma \left(\overline{r}\right)J_0\left(k\overline{r}\right)$$

(11)

and make use of the electrostatic potentials:

$${\mathsf{\Phi }}_1(r,z)=\underset{0}{\overset{\infty }{\int }}dkka\left(k\right)J_0\left(kr\right)e^{-\sqrt{k^2+\frac{1}{l_D^2}}z},{\mathsf{\Phi }}_2(r,z)=\underset{0}{\overset{\infty }{\int }}dkka\left(k\right)J_0\left(kr\right)e^{kz}.$$

(12)

Inserting Equations (11) and (12) into the boundary condition specified in Equation (4) allows us to calculate the zeroth order Hankel transform:

$$a\left(k\right)=\frac{1}{{\epsilon }_0}\frac{\underset{0}{\overset{\infty }{\int }}drr\sigma \left(r\right)J_0\left(kr\right)}{{\epsilon }_w\sqrt{k^2+\frac{1}{l_D^2}}+{\epsilon }_{l}k}=\frac{1}{{\epsilon }_0}\frac{\overline{\sigma }{r}_0{J}_1\left(k{r}_0\right)}{k\left({\epsilon }_w\sqrt{k^2+\frac{1}{l_D^2}}+{\epsilon }_{l}k\right)},$$

(13)

of the surface potential ${\mathsf{\Phi }}_s^{self}\left(r\right)={\mathsf{\Phi }}_1(r,z=0)={\mathsf{\Phi }}_2(r,z=0)$, where the equality on the right-hand side of Equation (13) is based on using the surface charge density in Equation (2) with $R\to \infty$. Knowing the function $a\left(k\right)$, introducing the dimensionalized wave number $y=kR$, and replacing $\overline{\sigma }={\sigma }_0{R}^2/r_0^2$ allow us to express the scaled surface potential as:

$$\frac{{\mathsf{\Phi }}_s^{self}\left(r\right)}{{\mathsf{\Phi }}_0}=\frac{R}{r_0}\underset{0}{\overset{\infty }{\int }}dy\frac{J_0\left(y\frac{r}{R}\right)J_1\left(y\frac{r_0}{R}\right)}{\sqrt{{\left(y\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}y}.$$

(14)

Note that ${\mathsf{\Phi }}_s\left(r\right)$ in Equation (10) applies to a circular disk of size $r_0$ and charge density $\overline{\sigma }$ inside a unit cell of radius R, whereas ${\mathsf{\Phi }}_s^{self}\left(r\right)$ in Equation (14) describes the surface potential in the limit on an infinitely large unit cell for the same circular disk of size $r_0$ and charge density $\overline{\sigma }$. In both cases, R, $r_0$, and ${\sigma }_0$—which is contained in ${\mathsf{\Phi }}_0={\sigma }_0{l}_D/\left({\epsilon }_0{\epsilon }_w\right)$—define the surface charge density $\overline{\sigma }={\sigma }_0{R}^2/r_0^2$ of the charged disk. The left diagram of Figure 4 shows an example for ${\mathsf{\Phi }}_s^{self}\left(r\right)/{\mathsf{\Phi }}_0$ (black curve), calculated for $r_0/l_D=1$, $R/l_D=1$, and ${\epsilon }_l={\epsilon }_w$. That is, the isolated charged disk described by the black curve has a radius equal to the Debye length and a surface charge density ${\sigma }_0$. Decreasing the unit cell radius from infinity to finite values (while keeping the disk unchanged) increases the scaled potential in the unit cell, ${\mathsf{\Phi }}_s\left(r\right)/{\mathsf{\Phi }}_0$, as shown for $R/l_D=4$ (purple), $R/l_D=2$ (red), $R/l_D=1.5$ (green), $R/l_D=1.2$ (orange), and $R/l_D=1$ (blue). In the final case (the blue line), the sizes of unit cell and disk are identical, $r_0=R=l_D$, implying ${\mathsf{\Phi }}_s\left(r\right)/{\mathsf{\Phi }}_0=1$.

The right diagram of Figure 4 shows the scaled difference potential $▵{\mathsf{\Phi }}_s\left(r\right)/{\mathsf{\Phi }}_0$, defined through $▵{\mathsf{\Phi }}_s\left(r\right)={\mathsf{\Phi }}_s\left(r\right)-{\mathsf{\Phi }}_s^{self}\left(r\right)$, for each displayed case in the left diagram, with the same color coding. For example, the blue curve in the right diagram displays the difference between the blue and black curves in the left diagram. Importantly, $▵{\mathsf{\Phi }}_s\left(r\right)$ corresponds to the electrostatic potential produced by all charged disks other than the central one located in our unit cell. Below, we will use $▵{\mathsf{\Phi }}_s\left(r\right)$ to calculate the Debye-Hückel free energy of the electric double layer at the dielectric interface.

We also point out that the limit $r_0\to 0$ of Equation (14) describes a single point charge that is located at the dielectric interface. The resulting scaled surface potential:

$${\left(\frac{{\mathsf{\Phi }}_s^{self}\left(r\right)}{{\mathsf{\Phi }}_0}\right)}_{r_0\to 0}=\frac{1}{2}\underset{0}{\overset{\infty }{\int }}dyy\frac{J_0\left(y\frac{r}{R}\right)}{\sqrt{{\left(y\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}y}$$

(15)

was first computed by Stillinger [22] and was later extended to off-interfacial locations of the charge and to bilayer geometries [24,41,42]. Hurd [43] demonstrated that in the limit ${\epsilon }_l\ll {\epsilon }_w$ and $l_D\ll r$, Equation (15) can be expressed as:

$${\left(\frac{{\mathsf{\Phi }}_s^{self}\left(r\right)}{{\mathsf{\Phi }}_0}\right)}_{r_0\to 0}=\frac{1}{2}\frac{R^2}{l_D^2}\left[\frac{l_D}{r}{e}^{-r/l_D}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D^3}{r^3}\right].$$

(16)

This is the sum of a screened Coulomb and dipole potential. The latter is long-range and is expected to dominate for sufficiently small average surface charge densities ${\sigma }_0=e/\left(\pi R^2\right)$.

2.3. Debye-Hückel Free Energy

Consider N isolated non-interacting circular disks, each of radius $r_0$ and located at the dielectric interface. Assembling the disks into an array of average surface charge density ${\sigma }_0$ and radius R of the unit cell is associated with an electrostatic interaction free energy per unit area $F/A$. We calculate that energy by multiplying the surface charge density at each point in a unit cell $\sigma \left(r\right)$ with the electrostatic potential produced by all other unit cells $▵{\mathsf{\Phi }}_s\left(r\right)$,

$$\frac{F}{A}=\frac{1}{2}\frac{\underset{0}{\overset{R}{\int }}drr\sigma \left(r\right)▵{\mathsf{\Phi }}_s\left(r\right)}{\underset{0}{\overset{R}{\int }}drr},$$

(17)

where the factor $1/2$ avoids double-counting of the cell-cell interactions. The free energy $F=F\left(r_0\right)$ in Equation (17) is a function of disk size, with the limit $r_0\to 0$ corresponding to discrete point-like charges. In the other limit, $r_0=R$, all charge is uniformly spread over the surface. In this case, $F(r_0=R)$ corresponds to the energy cost of merging initially isolated, preformed disks (all of radius R and surface charge density ${\sigma }_0$) into a uniformly charged interface with surface charge density ${\sigma }_0$. Note that $F(r_0=R)$ would become identical to $F_0$ upon ignoring the self-energy of the preformed disks (that is, upon replacing $▵{\mathsf{\Phi }}_s\left(r\right)$ by ${\mathsf{\Phi }}_s\left(r\right)$ in Equation (17)). Using the definition of $▵{\mathsf{\Phi }}_s\left(r\right)$, as well as the expressions for $\sigma \left(r\right)$ in Equation (2), for ${\mathsf{\Phi }}_s\left(r\right)/{\mathsf{\Phi }}_0$ in Equation (9), and for ${\mathsf{\Phi }}_s^{self}\left(r\right)/{\mathsf{\Phi }}_0$ in Equation (14), and integrating over the Bessel function give rise to the scaled Debye-Hückel free energy:

$$\frac{F}{F_0}=1+{\left(\frac{R}{r_0}\right)}^2\left\{\sum _{n=1}^{\infty }\frac{{\left[\frac{2J_1\left(y_n\frac{r_0}{R}\right)}{y_n{J}_0\left(y_n\right)}\right]}^2}{\sqrt{{\left(y_n\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}{y}_n}-\underset{0}{\overset{\infty }{\int }}dy\frac{\frac{2}{y}{\left[J_1\left(y\frac{r_0}{R}\right)\right]}^2}{\sqrt{{\left(y\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}y}\right\},$$

(18)

where $F_0$ is specified in Equation (1). Equation (18) is a major result of the present work; it is the relative change of the Debye-Hückel free energy of an ensemble of uniformly charged disks (each disk of radius $r_0$ in a Wigner–Seitz cell of radius R) as compared to the continuum limit, where due to the absence of charge discreteness, self-energies become irrelevant and where all charge is smeared out to a uniform surface charge density. This relative change is fully characterized by specifying the three ratios: $r_0/R$, $l_D/R$, and ${\epsilon }_l/{\epsilon }_w$. To obtain a numerical estimate of $F/F_0$, we carry out the summation and integration up to an upper limit only:

$$\frac{F}{F_0}=1+{\left(\frac{R}{r_0}\right)}^2\left\{\sum _{n=1}^{n_{max}}\frac{{\left[\frac{2J_1\left(y_n\frac{r_0}{R}\right)}{y_n{J}_0\left(y_n\right)}\right]}^2}{\sqrt{{\left(y_n\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}{y}_n}-\underset{0}{\overset{y_{n_{max}}}{\int }}dy\frac{\frac{2}{y}{\left[J_1\left(y\frac{r_0}{R}\right)\right]}^2}{\sqrt{{\left(y\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}y}\right\}.$$

(19)

Here, $y_{n_{max}}$ denotes the $n_{max}$’th zero of $J_1\left(y\right)$; see Figure 2. In Figure 5, we show $F/F_0$ as a function of $r_0/R$ for fixed $R/l_D=1$ and ${\epsilon }_l={\epsilon }_w$. Different curves correspond to different $n_{max}$.

Clearly, for growing $n_{max}$, the numerical approximation for $F/F_0$ in Equation (19) approaches a limiting value for every non-vanishing disk size $r_0>0$, yet not for the point charge limit $r_0\to 0$. To obtain an accurate estimate even for point charges, we point out that because the derivative of $J_1\left(r_{0}y\right)/r_0$ with respect to $r_0$ at the position $r_0=0$ vanishes, the derivative:

$${\left(\frac{d(F/F_0)}{d{r}_0}\right)}_{r_0=0}=0$$

(20)

also vanishes. Hence, we obtain an accurate numerical value of $F/F_0$ in the limit $r_0\to 0$ by calculating $F/F_0$ at a sufficiently small non-vanishing position $r_0$ and for sufficiently large $n_{max}$. More specifically, we choose $r_0\ll min(R,l_D)$ and $n_{max}\gg R/r_0$. The example shown in Figure 5 applies to $R=l_D$, for which $r_0\lesssim 0.1R$ and $n_{max}\gtrsim 3R/r_0$ are convenient choices. We will use that method instead of the formal limit of Equation (18),

$$\underset{r_0\to 0}{lim}\frac{F}{F_0}=1+\sum _{n=1}^{\infty }\frac{1}{\sqrt{{\left(y_n\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}{y}_n}\frac{1}{{\left[J_0\left(y_n\right)\right]}^2}-\underset{0}{\overset{\infty }{\int }}dy\frac{y/2}{\sqrt{{\left(y\frac{l_D}{R}\right)}^2+1}+\frac{{\epsilon }_l}{{\epsilon }_w}\frac{l_D}{R}y},$$

(21)

which contains two diverging contributions and thus is unsuitable to be computed numerically.

A this point, we have acquired the ability to numerically compute the Debye-Hückel free energy of a planar dielectric interface with an ensemble of discrete point charges (given by Equation (21), yet calculated using Equation (19) for sufficiently large $n_{max}$ and small $r_0/R$). The influence of considering discrete charges versus a smeared uniform charge density is directly reflected by the ratio $F/F_0$. This ratio depends only on the two parameters $l_D/R$ and ${\epsilon }_l/{\epsilon }_w$. The left diagram of Figure 6 shows how $F/F_0$ varies as a function of $l_D/R$.

The cell radius $R=\sqrt{e/\left({\sigma }_0\pi \right)}$ must be chosen such that one elementary charge $e={\sigma }_0\pi R^2$ is contained inside a single unit cell. Different curves correspond to different ${\epsilon }_l/{\epsilon }_w$. The ratio approaches $F/F_0=0$ for sufficiently small $l_D$ and $F/F_0=1$ for sufficiently large $l_D$. Indeed, while for very small $l_D$, the individual charges are screened and thus do not significantly interact with each other, the discreteness of the charge becomes irrelevant at large Debye lengths, implying that F and $F_0$ become identical. Note also that $F/F_0$ increases with growing ${\epsilon }_l$. Recall that F represents the interaction energy of point charges that are initially separated on the dielectric interface and are then brought together to form an ensemble of average surface charge density ${\sigma }_0$. When ${\epsilon }_l=0$, the charges are only able to interact through the salt-containing medium, which eliminates long-range interactions. For ${\epsilon }_l>0$, an electric field also forms in the salt-free medium. This field is long-range and results from the dipole that each surface charge forms with its associated diffuse counterion cloud; see Equation (16). These additional long-range interactions cause $F/F_0$ to increase with growing ${\epsilon }_l$.

When changing the surface charge density ${\sigma }_0$ at fixed Debye length $l_D$, it is convenient to plot the scaled free energy $f_s=2l_D^3{\epsilon }_w{\epsilon }_{0}F/\left(e^{2}A\right)$ as a function of the scaled surface charge density $s=l_D^2\times ({\sigma }_0/e)$. The right diagram of Figure 6 shows a double-logarithmic plot of $f_s$ as a function of s, for different ratios of ${\epsilon }_l/{\epsilon }_w$. For the continuum model (see Equation (1), where the discrete nature of the charges is neglected), $f_s=s^2$, as indicated by the uppermost curve in purple (which has a slope of two in the double-logarithmic plot). The other curves correspond to different ratios of ${\epsilon }_l/{\epsilon }_w$. The lowest one, for ${\epsilon }_l/{\epsilon }_w=0$ (blue color), differs substantially from the second lowest one, derived for ${\epsilon }_l/{\epsilon }_w=0.01$ (orange color). Indeed, the absence of long-range interactions for ${\epsilon }_l/{\epsilon }_w=0$ causes $f_s$ to decrease exponentially, as dictated by the screened Coulomb interaction. For non-vanishing ${\epsilon }_l/{\epsilon }_w$, the charges on the dielectric interface exhibit dipolar interactions, which leads to the scaled free energy $f_s\sim s^{5/2}$. This follows immediately from $F/A\sim (N/A)\times 1/r^3$, as well as $N/A\sim {\sigma }_0$ and $1/r^2\sim {\sigma }_0$ plus the definitions of $f_s$ and s. Hence, we expect all curves with ${\epsilon }_l/{\epsilon }_w>0$ in the right diagram of Figure 6 to adopt a slope of $2.5$ for a sufficiently small scaled average surface charge density s. This is exactly what we observe as indicated by the short black line segment, which has a slope of $2.5$. The change from $f_s\sim s^2$ to $f_s\sim s^{2.5}$ as a function of decreasing s is the most decisive feature of the Debye-Hückel model for discrete surface charges at a dielectric interface. Its accurate quantification is the main achievement of the present work.

3. Conclusions

Charges are often located at dielectric interfaces between a salt-containing aqueous solution and a salt-free medium of low dielectric constant. When the surface density of the charges becomes sufficiently small, their discrete nature cannot be ignored. We derived a free energy expression (see Equation (21)) that is valid within the linear Debye-Hückel framework for an ensemble of discrete charges at a planar dielectric interface. When the charge-to-charge distance on the interface is much smaller than the Debye length, we recover the well-known Debye-Hückel free energy expression for a continuous charge distribution, with the free energy $F\sim {\sigma }_0^2$ being proportional to the square of the surface charge density ${\sigma }_0$. When the charge-to-charge distance on the interface becomes much larger than the Debye length, interactions between dipoles formed by the interfacial charges with their counterions in solution become dominant, changing the free energy to a relationship $F\sim {\sigma }_0^{2.5}$. Biological cells operate at a Debye length of $l_D\approx 1$ nm. The nearest-neighbor distances between charges on biomolecular aggregates are often significantly larger than the Debye length. For example, the nearest-neighbor distance between the charges on a lipid membrane with 10% anionic lipids is about three times larger than the Debye length, necessitating accounting for the discreteness of the charge when using free energies to model membrane properties, such as domain formation, bending stiffness, spinodal decomposition, charge reversal [44,45], and renormalization [46], or differential capacitance [47,48]. Our present work may assist in this kind of modeling.

We emphasize that our model is valid within the Debye-Hückel limit, where correlations and ion size effects are neglected. Applying Poisson–Boltzmann theory or other non-linear approaches that account for finite ion sizes to discrete charges at a dielectric interface will likely not yield analytical expressions, but can be implemented numerically.