Continuum theory of electrostatic correlations at charged surfaces

TL;DR

This paper develops a continuum electrostatic correlation theory at charged surfaces, extending existing models to better predict ion behavior in complex electrolytes beyond the Poisson-Boltzmann approximation.

Contribution

It introduces a boundary condition and scaling laws within the Bazant-Storey-Kornyshev framework, improving predictions for correlated electrolytes at charged interfaces.

Findings

Accurately reproduces Monte Carlo simulation results across coupling regimes.

Provides a predictive continuum model for correlated electrolytes.

Extensible to complex electrolyte and ionic liquid models.

Abstract

The standard model for diffuse charge phenomena in colloid science, electrokinetics and biology is the Poisson-Boltzmann mean-field theory, which breaks down for multivalent ions and large surface charge densities due to electrostatic correlations. In this paper, we formulate a predictive continuum theory of correlated electrolytes based on two extensions of the Bazant-Storey-Kornyshev (BSK) framework: (i) a physical boundary condition enforcing continuity of the Maxwell stress at a charged interface, which upholds the Contact Theorem for dilute primitive-model electrolytes, and (ii) scaling relationships for the correlation length, for a one-component plasma at a charged plane and around a cylinder, as well as a dilute z:1 electrolyte screening a planar surface. In these cases, the theory accurately reproduces Monte Carlo simulation results from weak to strong coupling, and extensions are possible for more complex models of electrolytes and ionic liquids.

Equations

$$\begin{array}{rl}G=& {\int }_{V}d\text{r}\left\{g+\rho \varphi -\frac{ϵ}{2}\left[{\left(\nabla \varphi \right)}^2+{\ell }_c^2{\left({\nabla }^2\varphi \right)}^2+\dots \right]\right\}\\ & +{\oint }_{S}d{\text{r}}_{\text{s}}{q}_s\varphi .\end{array}$$

$$ϵ({\ell }_c^2{\nabla }^2-1){\nabla }^2\varphi =\rho .$$

$${\mu }_i={\mu }_i^{\theta }+kT\ln (c_i)+z_{i}e\varphi +{\mu }_i^{\mathrm{ex}}$$

$$-\mathbf{f}=\nabla P=kT\sum _i\nabla c_i+\rho \nabla \varphi +\sum _i{c}_i\nabla {\mu }_i^{\mathrm{ex}},$$

$$\nabla \cdot {\text{uptau}}_e=\rho \mathbf{E}=\nabla \cdot (\widehat{ϵ}\mathbf{E})\mathbf{E}$$

$$\begin{split}\mathbf{\uptau}_{e}=&\epsilon\mathbf{EE}-\frac{1}{2}\epsilon\mathbf{E}^{2}\,\mathbf{I}+\epsilon{\ell_{c}}^{2}\Big{[}\left(\mathbf{E}\cdot\nabla^{2}\mathbf{E}\right)\,\mathbf{I}-\mathbf{E}\left(\nabla^{2}\mathbf{E}\right)\\ &-\left(\nabla^{2}\mathbf{E}\right)\mathbf{E}+\frac{1}{2}\left(\nabla\cdot\mathbf{E}\right)^{2}\,\mathbf{I}\Big{]},\end{split}$$

$$\nabla \cdot {\text{uptau}}_e=\rho \mathbf{E}-\frac{1}{2}{\mathbf{E}}^2\nabla ϵ+\frac{1}{2}(\nabla \cdot \mathbf{E}{)}^2\nabla (ϵ{{\ell }_c}^2).$$

$${\text{uptau}}_e-{\text{uptau}}_{e,\mathrm{surf}}=0.$$

$$\hat{n}\cdot\ell_{c}\nabla^{3}\phi=\nabla^{2}\phi\Big{\rvert}_{S}$$

$$P=-\frac{q_{s}^{2}}{2\epsilon}+kT\sum_{i}c_{i}\Big{\rvert}_{S}=-\hat{n}\cdot\mathbf{\uptau}_{e}\cdot\hat{n}+kT\sum_{i}c_{i}\Big{\rvert}_{S}.$$

$$R_{\mathrm{hole}}={\left(\frac{1}{6}\right)}^{1/4}{\left(\frac{ze}{q_s}\right)}^{1/2}$$

$${\ell }_c=\alpha R_{\mathrm{hole}}$$

$${\delta }_c={\alpha }_2{\left(\frac{z^2{\ell }_B}{{\ell }_{\mathrm{GC}}}\right)}^{{\alpha }_3}{\left(\frac{z^2{\ell }_B}{{\lambda }_D}\right)}^{{\alpha }_4},$$

$${\ell }_c\sim {{\ell }_B}^{1/4}(q_s/e{)}^{-1/8}{C_{\mathrm{ref}}}^{-1/6},$$

$${\delta }_c=0.35{\left(\frac{z^2{\ell }_B}{{\ell }_{\mathrm{GC}}}\right)}^{-1/8}{\left(\frac{z^2{\ell }_B}{{\lambda }_D}\right)}^{2/3}.$$

$$\Xi =2\pi z^3{{\ell }_B}^2{q}_s/e,$$

$${\ell }_{\mathrm{GC}}=e(2\pi z{\ell }_B{q}_s{)}^{-1}$$

$${{\alpha }_0}^2\Xi {\tilde{\nabla }}^4\tilde{\varphi }-{\tilde{\nabla }}^2\tilde{\varphi }=2\tilde{\rho }=2e^{-\tilde{\varphi }}$$

$$\begin{array}{rl}& \hat{n}\cdot ({{\alpha }_0}^2\Xi {\tilde{\nabla }}^3\tilde{\varphi }-\tilde{\nabla }\tilde{\varphi })=-2\\ & \hat{n}\cdot {\alpha }_0\sqrt{\Xi }{\tilde{\nabla }}^3\tilde{\varphi }={\tilde{\nabla }}^2\tilde{\varphi }\end{array}$$

$${\lambda }_D=\sqrt{\frac{ϵkT}{(z^2+z)e^2{C}_{\mathrm{ref}}}},$$

$${\delta }_c^2{\tilde{\nabla }}^4\tilde{\varphi }-{\tilde{\nabla }}^2\tilde{\varphi }=\tilde{\rho }=\frac{z{e}^{-z\tilde{\varphi }}-z{e}^{\tilde{\varphi }}}{z^2+z}.$$

$$\begin{array}{rl}& \hat{n}\cdot ({\delta }_c^2{\tilde{\nabla }}^3\tilde{\varphi }-\tilde{\nabla }\tilde{\varphi })={\tilde{q}}_s\\ & \hat{n}\cdot {\delta }_c\widetilde{{\nabla }^3\varphi }={\tilde{\nabla }}^2\tilde{\varphi }.\end{array}$$

$$C_D=\frac{ϵ}{{\lambda }_D}\cosh \left(\frac{{\tilde{\varphi }}_D}{2}\right),$$

$$\hat{n}\cdot [{\text{uptau}}_{e,1}-{\text{uptau}}_{e,2}]=\hat{n}\cdot [{\text{uptau}}_{e,1}-{\text{uptau}}_{e,2}{]}_{{\ell }_c=0}.$$

Full text

1 Introduction

Electrostatic correlations can significantly affect the structure and thermodynamic properties of the electrical double layer 1, 2, resulting in qualitative differences from mean-field Poisson-Boltzmann (PB) theory, such as like-charge attraction 3, 4 or over-screening of surface charge. Critical applications in biology, colloids, separations, or electrochemistry rely on or operate in the regime where correlation effects are critical.

Numerous models have been proposed to capture electrostatic correlations, typically with a complicated mathematical structure. Outwaithe derived a modified Poisson Boltzmann models to account for the fluctuation potential of a single ion interacting with a charged wall 5. The hypernetted chain approximation and mean spherical approximation closure to the Ornstein-Zernike equation involve solving integral equations involving the direct correlation functions of bulk charged spheres 6, 7, 8, 9 for the equilibrium structure. Further work based on classical Density Functional Theory determines equilibrium properties based on the minimization of an integro-differential free energy functional. Kierlik and Rosinberg implemented a model (termed the bulk fluid model 10) which captures correlations based on an perturbative expansion of density with direct correlation functions as an input from the mean spherical approximation 11. Voudkavinova et al. analyzed the accuracy of two other related density functional theories (reference fluid density 12, functionalized mean spherical approximation 13, and bulk fluid 11, 14) in comparison to Monte Carlo simulations, finding the reference fluid density approach to be most accurate 10. These theories implement the accurate Fundamental Measure Theory functional developed by Rosenfeld to describe excluded volume effects 15, 16. The additional electrostatic interactions beyond mean field are included in the excess free energy separately, without any modification to the mean-field electrostatic part of the free energy.

While these approaches often produce accurate density profiles, they can be involved to implement to a broad range of applications, for different geometries or dynamic problems, especially compared to the classical PB theory. To our knowledge, they also have not yet been shown to recover the correlated behavior of the counterion-only limit for counterions of infinitesimal size. A simpler mathematical structure could help with the application and interpretation of electrostatic correlations to a wider class of problems in physics, including electrokinetics, colloidal interactions, and electrochemical transport/reactions.

Bazant, Storey and Kornyshev (BSK) proposed a continuum framework to account for the nonlocal dielectric permittivity of ionic liquids resulting from ion-ion correlations 17 with a simple mathematical structure, building on intermediate coupling approximations of Santangelo 18 and Hatlo and Lue 19 for the one-component plasma. The model captures correlations based on expansions in terms of electric field, rather than ion density, in the free energy functional which leads to a higher-order Poisson equation. In so doing, the electrostatic correlations are included self-consistently in the definition of the electrostatic potential whose gradient determines the electrostatic force on an ion in the diffuse layer. The BSK theory provides a simple framework to predict charge density oscillations and over-screening phenomena in a variety of electrokinetic , electrochemical, biophysical, and colloidal phenomena in electrolytes and ionic liquids. The equations require a similar level of complexity to solve compared to the PB theory, which allows them to be applied to a broad group of applications. The theory was used to describe electrosmotic 20 and electrophoretic mobility 21 reversals in multivalent and concentrated electrolytes, as well as electroconvective instabilities in ionic liquids 22. It was also applied to the dynamics 23, 24, 25, 26, 27 and electrosorption 28, 29, 30, 31 at electrochemical interfaces for ionic liquids and concentrated solvent-in-salt electrolytes, including storage 32 and transport 33 in nanoporous media. Electrostatic correlations have a profound effect on colloidal interactions4, 34, 35, where they can induce like-charge attraction in multivalent electrolytes, also predicted by the BSK theory. The activity and solvation energy of electrolytes at high concentration was also studied including the electrostatic correlation effect 36, 37, 38, 39, 40 as well as the extent of ion pairing in confinement41. Finally, the electrostatic correlations given by the BSK theory were important in describing the conduction through biological ion channels 42, 43, 44, 45, 46. Despite the numerous applications, fundamental questions remain about the proper boundary conditions and correlation length required to complete the BSK theory.

Here, we show that the appropriate boundary constraint for the higher-order Poisson equation is based on an interfacial stress balance. With corrected boundary conditions, the BSK theory becomes exact in the strong and weak coupling limits for the one-component plasma, and agrees with Monte Carlo (MC) simulations in intermediate coupling. We also suggest scaling relationships for the correlation length without steric constraints in one-component plasma and in multivalent electrolytes, for all the scenarios in Fig. 1. We show how the correlation length can have a simple physical interpretation based on the correlation hole size of counterions at a charged surface. Although generalizations are possible, we restrict the analysis to a restricted primitive electrolyte with hard, spherical ions of equal size in a constant permittivity, $\epsilon$, medium and smeared out surface charge density, $q_{s}$, and neglect all concentrated-solution effects, so as to isolate electrostatic correlations.

2 Theory

The BSK free energy functional is given by:

[TABLE]

Here, $g=(H-TS)/V$ is the enthalpy and entropy density, $\rho$ is the charge density, and $\phi$ is the electrostatic potential. For simplicity, the free energy is truncated after the first correlation contribution, although higher order terms can be considered 19. While the original authors performed a gradient expansion to arrive at Eq. 1 17, 20, the mathematical procedure is equivalent to modifying the interaction potential between ions from $U_{\alpha\beta}(r)=z_{\alpha}z_{\beta}\ell_{B}\mid r\mid^{-1}$ to $U_{\alpha\beta}(r)=z_{\alpha}z_{\beta}\ell_{B}\mid r\mid^{-1}(1-e^{-\mid r\mid/\ell_{c}})$ 25, 18. The modified interaction potential is solved in the mean-field limit. Thus the BSK theory is a phenomenological correction to PB within the mean-field approximation, by subtracting out interactions with smeared out charges within a correlation length $\ell_{c}$, which should scale as the size of the correlation hole. The modified Poisson equation results by finding the extremal of the functional ($\frac{\delta G}{\delta\phi}=0$):

[TABLE]

Eq. 2 is a statement of Maxwell’s equation, $\nabla\cdot\mathbf{D}=\rho$ where the displacement field is $\mathbf{D}=\hat{\epsilon}\mathbf{E}$ with a non-local permittivity operator $\hat{\epsilon}=\epsilon(\ell_{c}^{2}\nabla^{2}-1)$ applied on the electric field, $\mathbf{E}=-\nabla\phi$, in a medium of constant permittivity, $\epsilon$. PB theory is recovered when $\ell_{c}=0$. Note that the definition of electrostatic potential itself has changed by adding the higher order correlation terms, without violating Maxwell’s equation. In other words, the potential that determines the energy of an ion in the double layer must satisfy $\nabla\cdot\left(\hat{\epsilon}\mathbf{E}\right)=\rho$ rather than $\nabla\cdot\left(\epsilon\mathbf{E}\right)=\rho$. In this way, the electrostatic correlation contribution to the electrostatic energy for an ion is included self consistently within the electrostatic framework, rather than being added as additional corrections in the excess chemical potential. An advantage of this approach is that the electrostatic force per unit charge of an ion is captured directly with $\mathbf{E}$, meaning that the diffuse potential $\phi$ here could be measured experimentally at an electrode (if also accounting for the potential drop in the Stern layer).

The charge density at equilibrium, $\rho=\sum_{i}z_{i}ec_{i}$, will be determined by the constraint that the electrochemical potential for each ion is a constant at equilibrium. The electrochemical potential can be defined as the variation of the Gibbs free energy with respect to concentration 47, $\mu_{i}=\frac{\delta G}{\delta c_{i}}$, or

[TABLE]

where the first term is a reference value, the second term is the ideal entropy contribution, the third term is the electrostatic potential contribution, and the fourth term is the excess electrochemical potential.

The first open question in applying BSK theory is that of additional boundary conditions, beyond Maxwell’s equation $\hat{n}\cdot\mathbf{D}=q_{s}$. Presumably, the boundary condition must take care of the unaccounted short-range part of $U_{\alpha\beta}$. In the original BSK formulation and all subsequent works, the boundary condition of $\hat{n}\cdot\nabla^{3}\phi=0$ was applied, with the justification that the correlation effects should disappear at the interface 20, 25, 24, 33, 41, 48, 21, 34, 42, 26, 28, 27, 29, 49. The theory provided reasonable agreement to simulation and experimental results for ionic liquids and multivalent electrolytes. However, the boundary conditions have not yet been proved or validated systematically.

The second open question is the choice of correlation length, which has been arbitrarily set to the Bjerrum length for electrolytes 17, 21, and the ion diameter for ionic liquids 17. The theory is ultimately very sensitive to the choice of boundary conditions and correlation length. Here, we analyze the boundary condition in terms of a stress balance at the interface and then validate $\ell_{c}$ by comparison to MC simulations.

2.1 Interfacial balance

Applying the Gibbs-Duhem equation at constant temperature to the electrolyte and screened surface charges, following 50, and neglecting the external electrostatic work done on the system, gives $dP=\sum_{i}c_{i}d\mu_{i}$. Taking the gradient in three-dimensional space and applying the definition of the electrochemical potential:

[TABLE]

where $\mathbf{f}$ is the total thermodynamic force. The first and third terms on the RHS of Eq. 4 correspond to the gradient of osmotic pressure, $\nabla\Pi$. For an ideal solution, $\mu_{i}^{\mathrm{ex}}$=0. The gradient of the defined thermodynamic pressure is equivalent to the divergence of the total stress tensor of the electrolyte system, $\mathbf{f}=\nabla\cdot{\mathbf{\uptau}}$. The total stress tensor is composed of an osmotic pressure component, $\Pi$ and a Maxwell stress tensor, $\mathbf{\uptau}_{e}$, such that $\mathbf{\uptau}=-\Pi\,\mathrm{I}+\mathbf{\uptau}_{e}$. The component of interest in this analysis, $\mathbf{\uptau}_{e}$, can be defined by:

[TABLE]

in a constant $\epsilon$ medium. Plugging in for the charge density using the BSK Eq. 2 and performing integration by parts, one arrives at an expression for the Maxwell stress tensor for a fluid with a non-local permittivity,

[TABLE]

as derived in the Supporing Information. While the above equation was derived for constant $\epsilon$ and $\ell_{c}$, the expression is identical if these parameters vary. For varying $\epsilon$ or $\ell_{c}$, the Korteweg-Helmholtz force density must be included in the electrostatic stress 51,

[TABLE]

Within the distance of closest approach of the ions to the surface, correlations cannot affect the value of the Maxwell stress at the surface, $\mathbf{\uptau}_{e,\mathrm{surf}}$, generated by the surface charge density. The mechanical equilibrium problem therefore requires continuity in the electrostatic stress tensor evaluated in the solution and at the surface,

[TABLE]

At a uniformly-charged, flat surface without a dielectric jump, the Maxwell stress tensor is simply $\mathbf{\uptau}_{e,\mathrm{surf}}={q_{s}}^{2}/(2\epsilon)\hat{n}\hat{n}$, while the Maxwell stress tensor in the electrolyte is given by Eq. 6. Equating these two expressions, and substituting in $\hat{n}\cdot\mathbf{D}=q_{s}$ we arrive at a final boundary for a potential varying in one coordinate direction:

[TABLE]

applied at the distance of closest approach of the ion with the wall.

The method amounts to applying the Contact Theorem to the correlated electrolyte in the absence of correlations, shown below for $\mu_{i}^{\mathrm{ex}}=0$ at a flat electrode with constant charge density without a dielectric jump 52, 53:

[TABLE]

The Contact Theorem is a statement of mechanical equilibrium, where the repulsive osmotic pressure contribution is balanced by the electrostatic attraction from the Maxwell stress. Without the constraint from Eq. 8, the BSK theory does not obey this simple relationship which should be valid even for dilute electrolytes in the primitive model 54, 55, 56. The procedure of ensuring continuity in the Maxwell stress can be repeated for any higher order $\hat{\epsilon}$ by equating the $\mathbf{\uptau}_{e}$ at successive orders of derivatives. The condition in Eq. 8 is also applicable to any extended electrolyte mean-field theory with arbitrary models of concentrated solution activity and solvent polarizability, including interactions with a soft wall. The approach may even be extended to media with non-local dielectric constant $\hat{\epsilon}$ driven by solvent polarization 57, 58.

2.2 Correlation length scaling

At a highly charged surfaces, the charge-charge correlations are dominated by the mutual repulsion of counterions at the interface, as demonstrated in the schematic in Fig. 1(a). The size of a correlation hole of counterions forming a Wigner crystal is characterized by a length scale $R_{\mathrm{hole}}$:

[TABLE]

For this work, we assume that the correlation length scales as the size of a correlation hole, determined by the surface charge denisty:

[TABLE]

with one parameter $\alpha$, which is considered to be a constant. We will demonstrate that the fitted scaling with $\alpha=0.50$ works well from the limit of zero reservoir concentration (one-component plasma) to more concentrated electrolytes. At high concentration or at low surface charge densities, this scaling argument breaks down, and the other length scales might dominate in the correlation length. For example, if the surface charge tends to zero, then the charge-charge correlations will be dominated by the Bjerrum length and characteristic mean spacing between ions given by the bulk concentration. At very large charge densities and for large ion sizes, where $R_{\mathrm{hole}}$ becomes comparable to $a$, the ion diameter, the ion size can dominate in determining the charge-charge correlations due to over-crowding effects. Also, if other length scales are introduced, such as surface curvature, the correlation length between ions can also be affected, as demonstrated with the one-component plasma around a thin charged cylinder.

In the Supporting Information, the non-arbitrary scaling of the correlation length is investigated by comparison to Grand Canonical Monte Carlo simulations from ref. 59 for a $z:1$ electrolyte. Using the Buckingham-$\Pi$ theorem, we know that the dimensionless correlation length is related to a power law relationship, at least for small or large values of the dimensionless variables. Here we choose four length scales from which we can construct three dimensionless groups: $\ell_{c}$, $\ell_{\mathrm{GC}}$, $z^{2}\ell_{B}$, and $\lambda_{D}$, such that the powerlaw relationship can be expressed as:

[TABLE]

The correlation length scaling that arises from the fitting procedure to the MC data is given by:

[TABLE]

with a fitted scale of:

[TABLE]

Note that the fitted exponents are expressed in terms of fractions to emphasize their relationship to the intrinsic length scales in the system.

A detailed investigation of how correlations are affected as a function of surface charge density, ion valency, concentration, ion size, and surface curvature could motivate a more nuanced scaling of the correlation length, based on direct analysis of the charge-charge correlation function, including variations in the correlation length as a function of the distance from a charged surface. In this work, we isolate the electrostatic correlation effects for a dilute electrolyte at highly charged surfaces. For the purposes of the analysis in the main text, we will assume the correlation length to be $\ell_{c}=0.50R_{\mathrm{hole}}$ for all the scenarios investigated in Fig. 1. For the results with the fitted correlation length scaling in 15, one can refer to the Supporting Information.

2.3 One-component plasma

Considering a system of point-like counterions neutralizing a uniformly charged surface, the importance of correlations is governed by a coupling constant,

[TABLE]

which is a measure of the correlation hole size, $R_{\mathrm{hole}}$, compared to the characteristic ion distance from the surface, the Gouy-Chapman length,

[TABLE]

such that $\Xi\sim{R_{\mathrm{hole}}}^{2}/{\ell_{\mathrm{GC}}}^{2}$. Here, $z$ is the ion valency, $\ell_{B}$ is the Bjerrum length, and $e$ is an elementary charge. In the weak coupling limit ($\Xi<<1$), PB theory is valid. In the strong coupling limit ($\Xi>>1$), counterions interact with the electric potential created by the surface since ion-surface interactions dominate 54, 60, 55.

Now we consider applying the mechanical constraint, starting with the one-component plasma of infinitesimally small size with $\mu^{\mathrm{ex}}=0$. The one-component plasma consists of a single mobile ionic species which neutralizes the charge of a smeared out surface charge density. We can non-dimensionalize lengths with the Gouy-Chapman length, the potential by the thermal voltage for the counterion, $\phi_{0}=\frac{kT}{ze}$, and the charge density by $\rho_{0}=2\pi\ell_{B}q_{s}^{2}e^{-1}$. Here, we assume that the correlation length scales with the size of a correlation hole at the surface, $\delta_{c}=\ell_{c}/\ell_{\mathrm{GC}}=\alpha_{0}\sqrt{\Xi}$ Using $\thicksim$ to denote non-dimensionalized variables:

[TABLE]

with boundary conditions of

[TABLE]

at $\tilde{x}=0$, where $\alpha_{0}=1.36\alpha$ is an order one constant proportional to $\alpha$. Therefore, the importance of the higher order derivative is governed by the coupling parameter, $\Xi$.

The solution to these equations is compared to the results of MC simulations in Fig. 2 for a one-component plasma screening a plane of charge. The BSK theory reproduces the behavior of the one-component plasma from weak coupling, in intermediate coupling, and matches the strong coupling limit with $\alpha=0.50$. Furthermore, the BSK theory with the boundary condition of $\hat{n}\cdot\nabla^{3}\phi=0$, represented by the dashed-dotted lines in Fig. 2a, does not accurately represent the data at intermediate or strong coupling.

We can also consider the one-component plasma between two charged surfaces of equal charge density with the same sign, confining the counterions in a gap of dimensionless distance $\tilde{d}$, as shown in Fig. 3a. In Fig. 3b, the pressure is plotted as a function of separation distances between two charged surfaces with different coupling parameters, using Eq. 8 and using the same value for $\alpha$. The BSK theory again provides good agreement with the results of the MC simulations at all the coupling parameters.

Another critical question is the validity of Eq. 8 at a curved interface. The simplest model system to test the hypothesis is the one-component plasma surrounding a charged cylinder of radius $\xi=R_{\mathrm{cyl}}/\ell_{\mathrm{GC}}$ at infinite dilution, corresponding to a cylindrical cell with outer radius $R_{\mathrm{out}}\rightarrow\infty$. As shown in Fig. 4, the BSK equation reproduces the results of the weak and strong coupling limits correctly. However, similar to the strong coupling expansion of 62, the theory does not correctly describe the renormalization of charge arising from Manning condensation in the needle limit, where a fraction $f=1-1/\xi$ of the charge is “condensed" onto cylinder 63. The charge density must be multiplied by this fraction, $f$, in order to match the strong coupling expansion taking into account the charge renormalization in the needle limit 61, 64. The smaller the radius of curvature, the more likely that the configuration of correlated ions is influenced by curvature. In the “needle limit," where $\sqrt{\Xi}/\xi>>1$, ions are distributed in a nearly linear fashion along the cylindrical backbone with spacing scaling as $\tilde{a}\thicksim\Xi/\xi$, which may also be the relevant scaling for the correlation length in this regime rather than $\sqrt{\Xi}$. Supporting Figs. S1 and S2 show the results choosing $\delta_{c}=\Xi/\xi$, but the counterion condensation transition is still not captured.

3 Electrolytes

A more useful and relevant application of the BSK theory is to describe the distribution of charges in electrolytes and ionic liquids, as was originally proposed. Here, we focus on the dilute electrolyte limit, to isolate electrostatic correlations directly, without complications from overcrowding.

If the BSK equation for a z:1 electrolyte with salt concentration ${C_{\mathrm{ref}}}$ is non-dimensionalized with the thermal voltage $\tilde{\phi}=(e\phi)/(kT)$ and the Debye length

[TABLE]

$\tilde{\nabla}=\lambda_{D}\nabla$ and $\delta_{c}=\ell_{c}/\lambda_{D}$ the BSK equation becomes:

[TABLE]

The boundary conditions are similarly modified to:

[TABLE]

The agreement of the predicted charge density profiles from Eq. 21 with the GCMC data is very good, as exhibited in Fig. 4 for a 0.1 M 2:1 electrolyte. In the Supporting Information, the results are expanded to a more complete set of comparisons with simulations. It is seen with $\ell_{c}=0.50R_{\mathrm{hole}}$ or with $\ell_{c}$ determined by Eq. 15, that the BSK theory can correct the PB charge density profiles, including an overscreening transition. Larger errors from the BSK theory are incurred at large concentration, where the current assumption of $\mu^{\mathrm{ex}}_{i}=0$ breaks down.

One implication of the boundary condition is that the differential capacitance for $\ell_{c}=0$ is equivalent to the case of $\ell_{c}\neq 0$ if $\mu_{i}^{\mathrm{ex}}=0$. Therefore, the differential capacitance for the correlated, dilute electrolyte is given by the traditional Gouy-Chapman equation:

[TABLE]

in stark contrast to the original work in the limit of $\mu_{i}^{ex}=0$ 17, 20. It would be interesting to explore the implications of the boundary condition on electrokinetic reversals, electrochemical interfaces, biological channels, or colloidal phenomena 20, 25, 24, 33, 41, 48, 21, 34, 42, 26, 28, 27, 29. For example, the DLVO theory of colloidal interactions can be modified to include attractive correlation effects 4.

4 Further extensions of the theory

4.1 Charged Dielectric Interfaces:

Note that the Maxwell stress condition (Eq. 6) has only been stated without a dielectric jump. The stress condition may need further validation at a dielectric interface. A more general statement of matching the Maxwell stress with and without correlations might be given by a jump condition between the two media:

[TABLE]

For a uniformly-charged, flat interface without a dielectric jump, $\tau_{e,2}=0$, so the RHS of the above equation reduces to $\hat{n}\cdot[\mathbf{\uptau}_{e,1}-\mathbf{\uptau}_{e,2}]_{\ell_{c}=0}=\frac{{q_{s}}^{2}}{2\epsilon_{1}}\hat{n}$ .

4.2 Concentrated Electrolytes and Ionic Liquids:

The present analysis attempts to isolate the effect of ion correlations in a dilute electrolyte. Ion size effects, particularly for $a/\lambda_{D}>1$ will require further validation to properly account for correlations guided by ion size combined with electrostatics. A non-local free energy functional might be necessary to capture the size correlations in concentrated solutions 16, 65, in conjunction with electrostatic correlations. Short-range bulk correlations 66 are not captured in this theory. Furthermore, if surface charges are discrete rather than smeared out, the contact condition may change 67. For an arbitrary mixture of ions with different valency, the effective correlation length will depend upon correlations between each pair of species, although the correlations at high surface charge density will still be dominated by the most highly charged counterion.

5 Conclusions

The phenomenological BSK theory describes non-local, discrete correlation effects with a higher-order, local, continuum description of the free energy quite well. The remarkable agreement of the theory with the one-component plasma and primitive model electrolyte suggest that higher-order, continuum equations can properly account for correlation effects, as long as the appropriate constraints are imposed at boundaries. The formalism used here could be extended to the ionic liquid limit, although ion pairing, short-range non-electrostatic correlations, and “spin glass" ordering 68 might preclude a simple continuum description. Further modifications are needed to capture the long range screening exhibited in ionic liquids and concentrated electrolytes69, as well as density oscillations expected in overcrowded systems 16. Even so, the BSK theory captures important features of electrostatic correlations, including like-charge attraction and overscreening, driven by electrostatic interactions of spatially correlated counterions. Furthermore, unlike many previous approaches, all the electrostatic forces are contained self-consistently within the electrostatic potential, $\phi$. Detailed analysis of experimental data is needed to determine the competing effects of surface adsorption reactions modifying fixed surface charge 70 or the overscreening/ like-charge attraction effects modeled by the BSK theory 4.

6 Acknowledgements

This work was supported by an Amar G. Bose Research Grant (electrolyte calculations) and the Center for Enhanced Nanofluidic Transport (one-component plasma calculations). JPD is also supported by the National Science Foundation Graduate Research Fellowship under Grant No. 1122374. JPD would like to acknowledge useful discussions with Michael McEldrew, Amir Levy, Tingtao Zhou, Dimitrios Fraggedakis, and Mohamed Mirzadeh.