LETTER pubs.acs.org/JPCL
Free-Energy Density Functional of Ions at a Dielectric Interface Dirk Gillespie* Department of Molecular Biophysics and Physiology, Rush University Medical Center, Chicago, Illinois 60612, United States ABSTRACT: Ions at dielectric interfaces are found in a wide range of applications including biology, nanofluidics, and fuel cells. Often, the excluded volume of the ions has first-order effects on measured properties. Here, density functional theory of fluids is used to develop a statistical mechanical theory of how ions (charged hard spheres of different sizes) arrange around a dielectric interface. SECTION: Statistical Mechanics, Thermodynamics, Medium Effects
I
ons at a polarizable interface are common in many areas of science. Examples include ions (e.g., Naþ and Ca2þ in water or polyvalent polymers in ionic liquids) near the interface of two immiscible fluids, metal electrodes, proteins, lipid bilayers, or silica.15 While particle simulations (e.g., Monte Carlo, MC) can be used to compute the structure of the electrical double layer at a dielectric interface,4,6 there are few thermodynamic theories that directly compute the ensemble-averaged density profiles of the ions. While the interaction of two ions via a (linear) dielectric interface is well-understood, the resulting ion interaction is fundamentally different from the Coulomb potential because it is not longer radial. This makes the statistical mechanical averaging of particle interactions into ensemble-averaged quantities challenging. Recently, such a theory has been derived for point charges.7,8 In many nanoscale applications, however, the ions’ excluded volume has first-order effects.13 Outhwaite and Bhuiyan9 derived a modified PoissonBoltzmann equation for hard-sphere ions at a dielectric interface, but this is applicable only to the planar geometry and to ions of the same size. In this Letter, classic density functional theory (DFT) of charged hard-sphere ions with different sizes is used to derive a general theory of the dielectric electrical double layer that can be extended beyond the planar geometry tested here. Even in the planar geometry, it is the first thermodynamic theory of different-sized ions at a dielectric interface and the first DFT for dielectric interfaces (to my knowledge). In DFT of fluids, equilibrium density profiles {Fk(x)} are computed by minimizing an approximate free-energy functional, usually by solving the corresponding EulerLagrange equations.10 Accurate DFTs of ions in the primitive model of electrolytes (i.e., charged, hard spheres in a background dielectric) exist when the system has a single, uniform dielectric constant.1114 This kind of DFT is the starting point for the DFT with a dielectric interface. These Coulombic DFTs are generally based r 2011 American Chemical Society
on approximations of the direct correlation function and are derived either by functional integration11 or functional Taylor series.12,13 The results shown here use a generalization of the latter method14 that has been shown to be accurate compared to MC over a wide range of conditions.15 To build the dielectric functional, throughout this Letter, we consider ions in a region with dielectric constant ɛions near a dielectric interface S with dielectric constant ɛ. The ions cannot move across S and are in equilibrium with a bath electrolyte with concentrations {Fbath k }. For ions of species i, their radius is Ri and valence is zi. In a convention used throughout, suppose that an ion of species i is at x, and we consider its effect on a test ion of species j at x0 . At each s ∈ S, the presence of a charge zie at x (e is the fundamental charge) induces a surface charge σi given by6,16 σ i ðs; x, ^E Þ þ
Z ^E ^E zi e nðsÞ 3 ðs xÞ nðsÞ 3 ðs s0 Þ 0 σi ðs0 ; x, ^E Þ ds ¼ 3 0 2π S 2π Eions js xj3 js s j
ð1Þ where ^ɛ (ɛions ɛ)/(ɛions þ ɛ) and n(s) is the unit vector normal to S pointing out of the wall dielectric into the ion dielectric. The two ions then “feel” each other through the 2 C 0 0 potential ψCij (x,x0 ) þ ψD ij (x,x ), where ψij (x,x ) = zizje / 0 4πɛionsɛ0|x x | is the Coulomb potential with dielectric constant ɛions and Z zj e σ i ðs; xÞ ds ð2Þ ψDij ðx, x 0 Þ ¼ 4πE0 S js x 0 j (In some special geometries, both this reaction potential and the induced surface charge can be computed with image charges, but Received: February 8, 2011 Accepted: April 26, 2011 Published: April 28, 2011 1178
dx.doi.org/10.1021/jz2001908 | J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER
the theory presented here does not make use of image charges.) The ion’s effect on itself is included through its self-energy ψD ii (x, x), which depends only the ion’s position and can be treated as an external potential acting on the ion (see the second term in eq 16 later). 0 The potential ψD ij (x,x ) is not, however, part of the interaction potential between the ions; the ions themselves still interact only through the Coulomb potential (unlike, for example, if a Lennard-Jones interaction were added on top of the Coulomb). Instead, the dielectric reaction potential acts as an external potential, albeit an unusual one that depends on the ions’ densities, because it requires the creation of the induced surface charge σi(s;x,ɛ) on the geometric interface. Adding a dielectric interface is then more like adding a Lennard-Jones interaction to the hard wall. 0 If ψD ij (x,x ) were considered part of the interaction potential, then there is a thermodynamically exact way to construct a functional in DFT given by10 Z 1 ZZ 1 ð2Þ Fij ðx, x 0 ; RÞψDij ðx, x 0 Þ dx dx0 dR F D ½fFk ðxÞg ¼ 2 i, j 0 Z 1 Fi ðxÞψDii ðx, xÞ dx þ ð3Þ 2 i
∑
∑
0 In eq 3, F(2) ij (x,x ;R) is the two-body distribution function for an ion of species i at x and an ion of species j at x0 with the interaction 0 potential ψRij (x,x0 ) = ψCij (x,x0 ) þ RψD ij (x,x ). When considering a perturbation of the interaction potential, one wishes to find the intrinsic Helmholtz free energy, the part independent of the 0 external potential. In that case, the densities in F(2) ij (x,x ;R) would be {Fk(y)} of the R = 1 fully perturbed system. This is done by formally varying the external potential with R to keep the 17 0 densities in F(2) ij (x,x ;R) independent of R. In our case, however, we consider R as charging up the dielectric external potential; the same arguments that give eq 3 for interaction potentials18 give it for external potentials. Written in terms of the total correlation function (TCF), this gives ð2Þ
Fij ðx, x 0 ; RÞ ¼ Fi ðx; RÞFj ðx 0 ; RÞðhij ðx, x 0 ; RÞ þ 1Þ
ð4Þ
The density profiles {Fk(y;R)} now explicitly depend on R because we are now considering the total Helmholtz free energy, not just the intrinsic component. In this way, we account in part for a change in the intrinsic Helmholtz free energy but also for the R-coupling with the external potential that is generated by the density distribution. We assume that the densities Fk change linearly with R so that Æ0æ Fk ðx; RÞ ¼ ðFk ðxÞ FÆ0æ k ðxÞÞR þ Fk ðxÞ
ð5Þ
where Fk(x) is the density profile that we aim to calculate (R = 1) and FÆ0æ k (x) is the density profile of the unperturbed system in the absence of dielectric interfaces with R = 0 (which is precomputed before the dielectric calculation using the functional described in ref 14). That this is a good approximation was checked with the MC simulations of the results shown later. Approximating hij(x,x0 ;R) is the crux of the dielectric functional. To do that, we assume that hij is a combination of Additive, Rescaled-to-be-charge-neutral, Multiplied-by-the-wall-Boltzmannfactor Yukawa functions, each part of which will be described below. This approximation of the TCF gives rise to the name ARMY dielectric functional. We further assume that the TCF has
only three components, (1) the screening ion cloud around an 0 ion (hion ij (x,x )), (2) the ion cloud around the hard wall in the 0 absence of a dielectric interface (hwall ij (x,x )), and (3) the ion cloud around the charge induced on the dielectric interface by an 0 ion and its screening ion cloud (hD ij (x,x ;R)). The first two are independent of dielectrics, while the third depends on R. By the additivity assumption 0 wall 0 D 0 hij ðx, x0 ; RÞ hion ij ðx, x Þ þ hij ðx, x Þ þ hij ðx, x ; RÞ
ð6Þ
when |x x0 | > Rij Ri þ Rj and hij(x,x0 ;R) = 1 otherwise (i.e., when the ions overlap). The ion TCF hion ij describes how ions of species j pack around a central ion of species i at x. To approximate this, we take the TCF 0 hbulk ij (|x x |) of the homogeneous bulk electrolyte with concentrations {Fbath k }. However, this screening cloud can extend past the dielectric interface, where ions are excluded. To prevent is multiplied by the Boltzmann factor that defines the this, hbulk ij accessible volume of species j (produced by the external potential 0 Vext j (x )) 0 0 bulk 0 hion ij ðx, x Þ qi ðxÞχj ðx Þhij ðjx x jÞ
ð7Þ
0 if |x x0 | g Rij and 1 otherwise, where χj(x0 ) = exp(Vext j (x )/ kT). This ion cloud is no longer charge-neutral; therefore, it is rescaled to be with zi e Z qi ðxÞ ¼ 0 0 zj eFbath hion j ij ðjx x jÞ dx
∑j
jx x 0 j g R ij
where the denominator is the charge of all species j around an ion of species i located at x. The wall TCF hwall ij is the ion concentration profile when R = 0, which is 0 when |x x0 | < Rij and otherwise Æ0æ 0 0 bath 1 hwall ij ðx Þ ¼ Fj ðx Þ=Fj
ð8Þ
This term would include, for instance, the normal repulsion of ions at uncharged walls that occurs even in the absence of dielectric interfaces (see, e.g., the blue lines in Figure 2A) and the attraction/repulsion of ions by a fixed surface charge (a case not considered here to focus solely on dielectric effects). 0 Lastly, we approximate the dielectic TCF hD ij (x,x ;R). Because D hij depends on the parameter R, it is insightful to better understand the relationship between R and ɛ. Because R ramps up the contribution of the dielectric reaction potential ψD ij in the intermediate potential ψRij , intuitively, one might expect that each R corresponds to a different wall dielectric constant ɛR between the limits ɛ0 = ɛions (no dielectric contribution of the wall at R = 0) and ɛ1 = ɛ (full dielectric contribution at R = 1). Equation 1, an exact equation for σi(s;x,ɛ^) allows us to (approximately) derive such a relationship. Over the past decade, much work has been done on approximating this equation for fast applications (reviewed, for example, by Bardhan et al.16). This work shows that there are several good approximations for σi(s;x,ɛ^) that separate the dielectric variable ^ɛ from the geometry of the dielectric interface so that σi(s;x,ɛ^) ≈ ξ(ɛ^)σi*(s;x) for some functions ξ and σi* (e.g., the BIBEE/P approximation16 with ξ(ɛ^) = ^ɛ, a separation that is exact for a single planar interface with σ*i given in standard textbooks19). 0 ^) by eq 2. The dielectric With such a separation, ψD ij (x,x ) ξ(ɛ reaction potential contribution to the intermediate potential ψRij D is RψD ^). On the other hand, if this ij , and therefore, Rψij Rξ(ɛ 1179
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER
intermediate reaction potential is due to a dielectric constant ɛR, ^R) by eq 2 (where ^ɛR = (ɛions ɛR)/(ɛions þ then RψD ij ξ(ɛ ɛR)). Equating these proportionality factors defines the correspondence R ≈ ξ(ɛ^R)/ξ(ɛ^) for ɛR ranging between ɛ and ɛions; that is, there is a one-to-one relationship between R and intermediate dielectric constants ɛR. (Note that we have only used the existence of ξ to establish the existence of the ɛR; the exact form of ξ is not required for this derivation nor for the 0 results below that approximate how hD ij (x,x ;R) depends on R.) With this one-to-one relationship in mind, the dielectic TCF 0 hD ij (x,x ;R) is the ion cloud that accumulates at a wall with dielectric constant ɛR because of the surface charge induced by an ion of species i at x and the ions that screen it. To approximate the effective charge of this screened ion, we consider two limits. First, when x is far from the dielectric interface, the screening shell is complete, and the effective charge is 0. Second, when the ion is touching the dielectric interface, the ion’s charge is so close that it dominates the induced charge because the finite size of the ions keeps other ions at least a distance of Rij away; the screened ion’s effective charge is zi. The change in effective charge between these two limits occurs as the screening ions fill in the remaining void between the dielectric interface and the central ion, namely, the cylinder of radius Ri and height |sx x| (sx is the point on S closest to x). Here, we approximate the distance dependence as changing exponentially with the Debye length λ, namely, exp((Ri |sx x|)/λ). While an approximation (which is validated a posteriori in the figures), it is the kind of distance dependence that one would expect by integrating the counter charge from linearized PoissonBoltzmann (LPB) theory in that cylinder. The effective surface charge then is ER Þ σeff i ðs; x, ^
ðRi jsx xjÞ=λ
¼ σi ðs; x, ^ER Þe
ð9Þ
because σi is proportional to the ion’s effective charge (eq 1). Next, to examine how ions accumulate around such a surface charge, we consider the screened potential around each surface element s. Because each individual surface element is small, we approximate this potential as a Yukawa potential (like in LPB theory) so that γ ðxÞ σeff ER Þ js x0 j=λ i ðs; x, ^ e φi ðs, x ; x, ^ER Þ i 0 4πE0 js x j 0
ð10Þ
with an as-yet-undetermined coefficient γi(x). Because the charge on a single surface element is infinitesimal, the potentials from all of the elements are (approximately) additive. Therefore, the potential from the entire surface charge at x0 is Z Φi ðx 0 ; x, ^ER Þ φi ðs, x 0 ; x, ^ER Þ ds ð11Þ S
0
With another LPB approximation, for |x x | g Rij zj e ext 0 hDij ðx, x 0 ; RÞ Φi ðx0 ; x, ^ER ÞeVj ðx Þ=kT kT
ð12Þ
and 0 when ions overlap. As before, this TCF is multiplied by the Bolzmann factor of the ions’ accessible volume, and it is rescaled so that the ion cloud is charge-neutral, giving Z 4πEE0 kT σi ðs; x, ^E Þ ds S Z γi ðxÞ ¼ σ E Þ js x0 j=λ i ðs; x, ^ bath 2 e χj ðx0 Þ dx 0 e2 zj Fj 0j js x 0 jx x j g Rij j
Note that while the arrangement of ions around a screened surface charge requires the most approximations, we show their validity a posteriori by comparing the resulting DFT with MC simulations (see the figures later). 0 To understand how hD ij (x,x ;R) depends on R, we again use the existence of the separation σi(s;x,ɛ^R) ≈ ξ(ɛ^R) σi*(s;x). There0 ^R) (by eq 9), so that fore, hD ij (x,x ;R) ξ(ɛ hDij ðx, x 0 ; RÞ RhDij ðx, x 0 Þ
ð13Þ
D 0 0 where hD ij (x,x ) = hij (x,x ;1). (Note that the existence of the factoring of σi(s;x,ɛ^R) is used only to estimate the dependence of D 0 0 hD ij (x,x ;R) on R. In the actual calculation of hij (x,x ), one should use the exact σi(s;x,ɛ^) from eq 1.) The DFT is derived by substituting eqs 48, 12, and 13 into eq 3. Integrating over R gives the functional Z Fi ðxÞψDii ðx, xÞ dx 2F D ½fFk ðxÞg ZZ þ ΔFi ðxÞΔFj ðx 0 ÞHij ðx, x 0 ; 3ÞψDij ðx, x 0 Þ dx dx 0 ZZ 0 0 D 0 0 þ ΔFi ðxÞFÆ0æ j ðx ÞHij ðx, x ; 2Þψij ðx, x Þ dx dx ZZ 0 D 0 0 þ ΔFj ðx0 ÞFÆ0æ i ðxÞHij ðx, x ; 2Þψij ðx, x Þ dx dx ZZ Æ0æ 0 0 D 0 0 þ FÆ0æ ð14Þ i ðxÞFj ðx ÞHij ðx, x ; 1Þψij ðx, x Þ dx dx
where, for space considerations, the Einstein summation convention is used, ΔFk(y) = Fk(y) FÆ0æ k (y), and 1 1 D 0 wall 0 h ðx, x 0 Þ Hij ðx, x0 ; nÞ ¼ ðhion ij ðx, x Þ þ hij ðx Þ þ 1Þ þ n n þ 1 ij ð15Þ For the functional derivative 1 Æ0æ δF D =δFi ðxÞ zi eðφÆ0æ CD ðxÞ φC ðxÞÞ 6 1 1 þ ψDii ðx, xÞ þ zi eðφCD ðxÞ φC ðxÞÞ 2Z 3 1 0 D 0 0 cj ðx 0 ; 2Þhion þ ij ðx, x Þψij ðx, x Þ dx 3 Z 1 0 D 0 0 cj ðx 0 ; 2Þhwall þ ij ðx, x Þψij ðx, x Þ dx 3 Z 1 cj ðx 0 ; 3ÞhDij ðx, x 0 ÞψDij ðx, x 0 Þ dx0 þ 4
ð16Þ
where 1 cj ðx 0 ; nÞ ¼ Fj ðx0 Þ þ FÆ0æ ðx 0 Þ n j
ð17Þ
0 This used the fact that, for ψij(x,x0 ) = ψCij (x,x0 ) þ ψD ij (x,x ) Z Fj ðx 0 Þψij ðx, x 0 Þ dx 0 ¼ zi eφCD ðxÞ ð18Þ
∑j
∑
∑j 1180
Z
Fj ðx 0 ÞψCij ðx, x 0 Þ dx0 ¼ zi eφC ðxÞ
ð19Þ
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER 19
where φC and φCD satisfy the Poisson equations E0 r 3 ½EðxÞrφCD ðxÞ ¼ e E0 Er2 φC ðxÞ ¼ e
∑j zj Fj ðxÞ
∑j zj Fj ðxÞ
ð20Þ ð21Þ
Æ0æ Æ0æ Æ0æ Similar definitions hold for φC and φCD using Fj (x0 ). Next, we compare the DFT calculations to MC simulations. In general, the DFT took ∼1 min on a single Macintosh PowerPC G5 2 GHz processor, which is substantially faster than the ∼1 h for the MC but slower than the Coulombic DFT, which uses fast Fourier transforms to compute convolutions created by the purely radial interaction potential.20 Figures 13 show the comparison of the ARMY dielectric functional at a single planar dielectric interface to MC simulations (done as previously described6). To focus solely on the effect of the dielectric interface, there is no additional fixed surface charge on the wall.
Figure 1. Concentration profiles of ions near a single planar dielectric interface at x = 0; ɛ is 2 (black), 20 (red), 40 (green), 80 (blue), and ¥ (magenta). Monovalent ions (diameter 3 Å) at bath concentrations of 0.01, 0.1, and 1 M. The dashed line shows the results of using only the self-energy term for ɛ = 2. In this and all figures, the lines are the DFT results, and the symbols are from MC simulations; the surface charge is 0, and T = 298.15 K, ɛions = 80, and ɛ = 500 were used in the DFT to approximate ¥.
Even this seemingly simple system is challenging because dielectric interactions are dominated by correlations. This can be seen by considering the restricted primitive model (RPM, where all ions have the same size) of monovalent ions. In this case, both cations and anions feel the same potential ψD ij because ψD ij zizj, and therefore, the cation and anion concentration profiles are identical and charge-neutral everywhere (Figure 1). If the correlations among ions are ignored in a mean-field approach, the mean-field averaged potential is identically 0 because the mean electrostatic potential φCD(x) = 0 everywhere, and therefore, δFD/δFi(x) = ψD ii (x,x)/2. This approximation is very poor, especially at high concentrations, as shown for ɛ = 2 with dotted lines in Figure 1. In this geometry, 1 19 and it is the hion ψD ii (x,x) decays as x , ij term in eq 16 that compensates for this far from the wall (not shown). This makes sense intuitively because this correlation describes the ion cloud around the central ion that produced the self-energy term. The hD ij term in eq 16 is important near the interface, sometimes becoming as large as the hion ij term (not shown). The contributerm in eq 16 was relatively small in all of the tion of the hwall ij cases studied (not shown) but will be larger when a fixed surface charge is included. The relative sizes of these terms then give some physical insight into the repulsion of ions by low dielectrics; the asympotic behavior is governed by the canceling of the self-energy and bulk screening cloud terms, while near the dielectric interface (∼8 Å for Figure 2A), the ion cloud near the induced surface charge and the distorted screening cloud around each ion dominate the formation of the double layer. The TCFs described above are based on principles derived from LPB, including additivity and eq 12. Therefore, it must be tested whether the functional fails when these approximations are questionable, namely, high valence and high concentration. This is done in Figures 2 and 3 for ion sizes and concentrations of normal electrolytes. Further testing of other conditions (e.g., significantly higher concentrations of ionic liquids5) will be reported in future work. Figure 2 shows the RPM for divalent cations and monovalent anions. Over the range of ɛ = 2¥ for = 0.5Fbath = 1 M, the ARMY dielectric functional works Fbath 2 1 extremely well. In fact, it reproduces many subtle features of the concentration profile found in the MC simulations. These include the steep initial slope (near x j 3 Å) in the anion profile, which then decreases, something not seen in the cation profile.
Figure 2. (A) Concentration profiles of divalent cations and monovalent anions (diameter 3 Å) at a 1 M bath cation concentration for ɛ = 2 (black), 20 (red), 40 (green), 80 (blue), and ¥ (magenta). (B) ɛ = ¥ for the same ions at 0.01 (black), 0.1 (blue), and 1 M (magenta) bath cation concentrations. Solid lines/shaded symbols: cations. Dashed lines/open symbols: anions. 1181
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER
’ ACKNOWLEDGMENT This material is based on work supported by, or in part by, the U.S. Army Research Laboratory and the U.S. Army Research Office under Contract No. W911NF-09-1-0488. Special thanks to Dezs€o Boda for providing the Monte Carlo simulations. Very useful discussions with Roland Roth and Jaydeep Bardhan are gratefully acknowledged. ’ REFERENCES
Figure 3. Concentration profiles of Kþ (solid lines/shaded symbols, diameter 2.76 Å), Ca2þ (dashed lines/open symbols, diameter 2 Å), and Cl (dotted lines/half-filled symbols, diameter 3.62 Å) at bath concentrations of 1, 0.1, and 1.2 M, respectively. The main part shows ɛ = ¥, and the inset shows ɛ = 2.
This is another physical insight: the monovalent ions actually have a shaper depletion near the dielectric interface than the divalents. The DFT also reproduces the small “bump” in the cation profile, except that the maximum of this bump is shifted to larger x. This and the tendencies for the DFT to be above the MC near the interface and below thereafter are not, however, defects of the dielectric functional. They mirror the same features in the ion profiles when there is no dielectric difference between the bath and the wall (ɛ = 80, blue lines). This is true for the vast majority of the profiles studied so far, indicating that many of the errors seen in the curves with dielectrics are due to the underlying Coulomb electrostatics functional. Figure 3 shows the results for ions of different size, specifically a mixture of Kþ, Ca2þ, and Cl . The main part of Figure 3 shows the results for ɛ = ¥. While the DFT tends to overestimate the adsorption of all ion species (a pattern also seen in Figure 1), the agreement with MC is generally good, with DFT reproducing the minimum in the Kþ profile and Ca2þ being adsorbed more than Kþ. The inset shows very good agreement for ɛ = 2. The “kink” in the Cl curve for ɛ = ¥ is an artifact for the largest ion only, stemming from the multiplication of bulk TCF by the Boltzmann factor of the wall geometry (data not shown). Its effect is small in the cases considered so far, but it is an area of improvement for the ARMY approach. Dielectric interfaces have first-order effects on the structure of ions that cannot be reproduced by mean-field approaches. The DFT derived here includes the necessary correlations to better understand the dielectric electrical double layer over a wide range of dielectric constants, ion concentrations, and ion species (including mixtures). In this first report, the planar geometry is considered, but the theory is general and will be tested in other geometries in future work. Even in the planar geometry, it is the first theory of ions with different sizes at a dielectric interface. Moreover, this work shows that a dielectric theory can be constructed as a perturbation of the Coulomb potential.
’ AUTHOR INFORMATION Corresponding Author
*E-mail: [email protected].
(1) Boda, D.; Fawcett, W. R.; Henderson, D.; Sokolowski, S. Monte Carlo, Density Functional Theory, and PoissonBoltzmann Theory Study of the Structure of an Electrolyte Near an Electrode. J. Chem. Phys. 2002, 116, 7170–7176. (2) \Boda, D.; Valisko, M.; Eisenberg, B.; Nonner, W.; Henderson, D.; Gillespie, D. Combined Effect of Pore Radius and Protein Dielectric Coefficient on the Selectivity of a Calcium Channel. Phys. Rev. Lett. 2007, 98, 168102. (3) van der Heyden, F. H. J.; Stein, D.; Besteman, K.; Lemay, S. G.; Dekker, C. Charge Inversion at High Ionic Strength Studied by Streaming Currents. Phys. Rev. Lett. 2006, 96, 224502. (4) Torrie, G. M.; Valleau, J. P. Double Layer Structure at the Interface Between Two Immiscible Electrolyte Solutions. J. Electroanal. Chem. 1986, 206, 69–79. (5) Loth, M. S.; Skinner, B.; Shklovskii, B. I. Anomalously Large Capacitance of an Ionic Liquid Described by the Restricted Primitive Model. Phys. Rev. E 2010, 82, 056102. (6) Boda, D.; Gillespie, D.; Nonner, W.; Henderson, D.; Eisenberg, B. Computing Induced Charges in Inhomogeneous Dielectric Media: Application in a Monte Carlo Simulation of Complex Ionic Systems. Phys. Rev. E 2004, 69, 046702. (7) Buyukdagli, S.; Manghi, M.; Palmeri, J. Variational Approach for Electrolyte Solutions: From Dielectric Interfaces to Charged Nanopores. Phys. Rev. E 2010, 81, 041601. (8) Hatlo, M. M.; Lue, L. Electrostatic Interactions of Charged Bodies from the Weak- to the Strong-Coupling Regime. Europhys. Lett. 2010, 89, 25002. (9) Outhwaite, C. W.; Bhuiyan, L. B. An Improved Modified PoissonBoltzmann Equation in Electric-Double-Layer Theory. J. Chem. Soc., Faraday Trans. 2 1983, 79, 707–718. (10) Evans, R. Density Functionals in the Theory of Nonuniform Fluids. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; pp 85176. (11) Tang, Z.; Scriven, L. E.; Davis, H. T. Interactions Between Primitive Electrical Double Layers. J. Chem. Phys. 1992, 97, 9258–9266. (12) Kierlik, E.; Rosinberg, M. L. Density-Functional Theory for Inhomogeneous Fluids: Adsorption of Binary Mixtures. Phys. Rev. A 1991, 44, 5025–5037. (13) Rosenfeld, Y. Free Enery Model for Inhomogeneous Fluid Mixtures: Yukawa-Charged Hard Spheres, General Interactions, and Plasmas. J. Chem. Phys. 1993, 98, 8126–8148. (14) Gillespie, D.; Nonner, W.; Eisenberg, R. S. Coupling Poisson NernstPlanck and Density Functional Theory to Calculate Ion Flux. J. Phys.: Condens. Matter 2002, 14, 12129. (15) Gillespie, D.; Valisko, M.; Boda, D. Density Functional Theory of the Electrical Double Layer: The RFD Functional. J. Phys.: Condens. Matter 2005, 17, 6609–6626. (16) Bardhan, J. P.; Knepley, M. G.; Anitescu, M. Bounding the Electrostatic Free Energies Associated with Linear Continuum Models of Molecular Solvation. J. Chem. Phys. 2009, 130, 104108. (17) Sullivan, D. E. Statistical Mechanics of a Nonuniform Fluid with Long-Range Attractions. Phys. Rev. A 1982, 25, 1669. (18) Davis, H. T. Statistical Mechanics of Phases, Interfaces, and Thin Films; Wiley-VCH: New York, 1996. (19) Jackson, J. D. Classical Electrodynamics, 3rd ed.; John Wiley and Sons, Inc.: New York, 1999. (20) Knepley, M. G.; Karpeev, D.; Davidovits, S.; Eisenberg, R. S.; Gillespie, D. An Efficient Algorithm for Classical Density Functional Theory in Three Dimensions: Ionic Solutions. J. Chem. Phys. 2010, 132, 124101. 1182
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
Free-Energy Density Functional of Ions at a Dielectric Interface Dirk Gillespie* Department of Molecular Biophysics and Physiology, Rush University Medical Center, Chicago, Illinois 60612, United States ABSTRACT: Ions at dielectric interfaces are found in a wide range of applications including biology, nanofluidics, and fuel cells. Often, the excluded volume of the ions has first-order effects on measured properties. Here, density functional theory of fluids is used to develop a statistical mechanical theory of how ions (charged hard spheres of different sizes) arrange around a dielectric interface. SECTION: Statistical Mechanics, Thermodynamics, Medium Effects
I
ons at a polarizable interface are common in many areas of science. Examples include ions (e.g., Naþ and Ca2þ in water or polyvalent polymers in ionic liquids) near the interface of two immiscible fluids, metal electrodes, proteins, lipid bilayers, or silica.15 While particle simulations (e.g., Monte Carlo, MC) can be used to compute the structure of the electrical double layer at a dielectric interface,4,6 there are few thermodynamic theories that directly compute the ensemble-averaged density profiles of the ions. While the interaction of two ions via a (linear) dielectric interface is well-understood, the resulting ion interaction is fundamentally different from the Coulomb potential because it is not longer radial. This makes the statistical mechanical averaging of particle interactions into ensemble-averaged quantities challenging. Recently, such a theory has been derived for point charges.7,8 In many nanoscale applications, however, the ions’ excluded volume has first-order effects.13 Outhwaite and Bhuiyan9 derived a modified PoissonBoltzmann equation for hard-sphere ions at a dielectric interface, but this is applicable only to the planar geometry and to ions of the same size. In this Letter, classic density functional theory (DFT) of charged hard-sphere ions with different sizes is used to derive a general theory of the dielectric electrical double layer that can be extended beyond the planar geometry tested here. Even in the planar geometry, it is the first thermodynamic theory of different-sized ions at a dielectric interface and the first DFT for dielectric interfaces (to my knowledge). In DFT of fluids, equilibrium density profiles {Fk(x)} are computed by minimizing an approximate free-energy functional, usually by solving the corresponding EulerLagrange equations.10 Accurate DFTs of ions in the primitive model of electrolytes (i.e., charged, hard spheres in a background dielectric) exist when the system has a single, uniform dielectric constant.1114 This kind of DFT is the starting point for the DFT with a dielectric interface. These Coulombic DFTs are generally based r 2011 American Chemical Society
on approximations of the direct correlation function and are derived either by functional integration11 or functional Taylor series.12,13 The results shown here use a generalization of the latter method14 that has been shown to be accurate compared to MC over a wide range of conditions.15 To build the dielectric functional, throughout this Letter, we consider ions in a region with dielectric constant ɛions near a dielectric interface S with dielectric constant ɛ. The ions cannot move across S and are in equilibrium with a bath electrolyte with concentrations {Fbath k }. For ions of species i, their radius is Ri and valence is zi. In a convention used throughout, suppose that an ion of species i is at x, and we consider its effect on a test ion of species j at x0 . At each s ∈ S, the presence of a charge zie at x (e is the fundamental charge) induces a surface charge σi given by6,16 σ i ðs; x, ^E Þ þ
Z ^E ^E zi e nðsÞ 3 ðs xÞ nðsÞ 3 ðs s0 Þ 0 σi ðs0 ; x, ^E Þ ds ¼ 3 0 2π S 2π Eions js xj3 js s j
ð1Þ where ^ɛ (ɛions ɛ)/(ɛions þ ɛ) and n(s) is the unit vector normal to S pointing out of the wall dielectric into the ion dielectric. The two ions then “feel” each other through the 2 C 0 0 potential ψCij (x,x0 ) þ ψD ij (x,x ), where ψij (x,x ) = zizje / 0 4πɛionsɛ0|x x | is the Coulomb potential with dielectric constant ɛions and Z zj e σ i ðs; xÞ ds ð2Þ ψDij ðx, x 0 Þ ¼ 4πE0 S js x 0 j (In some special geometries, both this reaction potential and the induced surface charge can be computed with image charges, but Received: February 8, 2011 Accepted: April 26, 2011 Published: April 28, 2011 1178
dx.doi.org/10.1021/jz2001908 | J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER
the theory presented here does not make use of image charges.) The ion’s effect on itself is included through its self-energy ψD ii (x, x), which depends only the ion’s position and can be treated as an external potential acting on the ion (see the second term in eq 16 later). 0 The potential ψD ij (x,x ) is not, however, part of the interaction potential between the ions; the ions themselves still interact only through the Coulomb potential (unlike, for example, if a Lennard-Jones interaction were added on top of the Coulomb). Instead, the dielectric reaction potential acts as an external potential, albeit an unusual one that depends on the ions’ densities, because it requires the creation of the induced surface charge σi(s;x,ɛ) on the geometric interface. Adding a dielectric interface is then more like adding a Lennard-Jones interaction to the hard wall. 0 If ψD ij (x,x ) were considered part of the interaction potential, then there is a thermodynamically exact way to construct a functional in DFT given by10 Z 1 ZZ 1 ð2Þ Fij ðx, x 0 ; RÞψDij ðx, x 0 Þ dx dx0 dR F D ½fFk ðxÞg ¼ 2 i, j 0 Z 1 Fi ðxÞψDii ðx, xÞ dx þ ð3Þ 2 i
∑
∑
0 In eq 3, F(2) ij (x,x ;R) is the two-body distribution function for an ion of species i at x and an ion of species j at x0 with the interaction 0 potential ψRij (x,x0 ) = ψCij (x,x0 ) þ RψD ij (x,x ). When considering a perturbation of the interaction potential, one wishes to find the intrinsic Helmholtz free energy, the part independent of the 0 external potential. In that case, the densities in F(2) ij (x,x ;R) would be {Fk(y)} of the R = 1 fully perturbed system. This is done by formally varying the external potential with R to keep the 17 0 densities in F(2) ij (x,x ;R) independent of R. In our case, however, we consider R as charging up the dielectric external potential; the same arguments that give eq 3 for interaction potentials18 give it for external potentials. Written in terms of the total correlation function (TCF), this gives ð2Þ
Fij ðx, x 0 ; RÞ ¼ Fi ðx; RÞFj ðx 0 ; RÞðhij ðx, x 0 ; RÞ þ 1Þ
ð4Þ
The density profiles {Fk(y;R)} now explicitly depend on R because we are now considering the total Helmholtz free energy, not just the intrinsic component. In this way, we account in part for a change in the intrinsic Helmholtz free energy but also for the R-coupling with the external potential that is generated by the density distribution. We assume that the densities Fk change linearly with R so that Æ0æ Fk ðx; RÞ ¼ ðFk ðxÞ FÆ0æ k ðxÞÞR þ Fk ðxÞ
ð5Þ
where Fk(x) is the density profile that we aim to calculate (R = 1) and FÆ0æ k (x) is the density profile of the unperturbed system in the absence of dielectric interfaces with R = 0 (which is precomputed before the dielectric calculation using the functional described in ref 14). That this is a good approximation was checked with the MC simulations of the results shown later. Approximating hij(x,x0 ;R) is the crux of the dielectric functional. To do that, we assume that hij is a combination of Additive, Rescaled-to-be-charge-neutral, Multiplied-by-the-wall-Boltzmannfactor Yukawa functions, each part of which will be described below. This approximation of the TCF gives rise to the name ARMY dielectric functional. We further assume that the TCF has
only three components, (1) the screening ion cloud around an 0 ion (hion ij (x,x )), (2) the ion cloud around the hard wall in the 0 absence of a dielectric interface (hwall ij (x,x )), and (3) the ion cloud around the charge induced on the dielectric interface by an 0 ion and its screening ion cloud (hD ij (x,x ;R)). The first two are independent of dielectrics, while the third depends on R. By the additivity assumption 0 wall 0 D 0 hij ðx, x0 ; RÞ hion ij ðx, x Þ þ hij ðx, x Þ þ hij ðx, x ; RÞ
ð6Þ
when |x x0 | > Rij Ri þ Rj and hij(x,x0 ;R) = 1 otherwise (i.e., when the ions overlap). The ion TCF hion ij describes how ions of species j pack around a central ion of species i at x. To approximate this, we take the TCF 0 hbulk ij (|x x |) of the homogeneous bulk electrolyte with concentrations {Fbath k }. However, this screening cloud can extend past the dielectric interface, where ions are excluded. To prevent is multiplied by the Boltzmann factor that defines the this, hbulk ij accessible volume of species j (produced by the external potential 0 Vext j (x )) 0 0 bulk 0 hion ij ðx, x Þ qi ðxÞχj ðx Þhij ðjx x jÞ
ð7Þ
0 if |x x0 | g Rij and 1 otherwise, where χj(x0 ) = exp(Vext j (x )/ kT). This ion cloud is no longer charge-neutral; therefore, it is rescaled to be with zi e Z qi ðxÞ ¼ 0 0 zj eFbath hion j ij ðjx x jÞ dx
∑j
jx x 0 j g R ij
where the denominator is the charge of all species j around an ion of species i located at x. The wall TCF hwall ij is the ion concentration profile when R = 0, which is 0 when |x x0 | < Rij and otherwise Æ0æ 0 0 bath 1 hwall ij ðx Þ ¼ Fj ðx Þ=Fj
ð8Þ
This term would include, for instance, the normal repulsion of ions at uncharged walls that occurs even in the absence of dielectric interfaces (see, e.g., the blue lines in Figure 2A) and the attraction/repulsion of ions by a fixed surface charge (a case not considered here to focus solely on dielectric effects). 0 Lastly, we approximate the dielectic TCF hD ij (x,x ;R). Because D hij depends on the parameter R, it is insightful to better understand the relationship between R and ɛ. Because R ramps up the contribution of the dielectric reaction potential ψD ij in the intermediate potential ψRij , intuitively, one might expect that each R corresponds to a different wall dielectric constant ɛR between the limits ɛ0 = ɛions (no dielectric contribution of the wall at R = 0) and ɛ1 = ɛ (full dielectric contribution at R = 1). Equation 1, an exact equation for σi(s;x,ɛ^) allows us to (approximately) derive such a relationship. Over the past decade, much work has been done on approximating this equation for fast applications (reviewed, for example, by Bardhan et al.16). This work shows that there are several good approximations for σi(s;x,ɛ^) that separate the dielectric variable ^ɛ from the geometry of the dielectric interface so that σi(s;x,ɛ^) ≈ ξ(ɛ^)σi*(s;x) for some functions ξ and σi* (e.g., the BIBEE/P approximation16 with ξ(ɛ^) = ^ɛ, a separation that is exact for a single planar interface with σ*i given in standard textbooks19). 0 ^) by eq 2. The dielectric With such a separation, ψD ij (x,x ) ξ(ɛ reaction potential contribution to the intermediate potential ψRij D is RψD ^). On the other hand, if this ij , and therefore, Rψij Rξ(ɛ 1179
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER
intermediate reaction potential is due to a dielectric constant ɛR, ^R) by eq 2 (where ^ɛR = (ɛions ɛR)/(ɛions þ then RψD ij ξ(ɛ ɛR)). Equating these proportionality factors defines the correspondence R ≈ ξ(ɛ^R)/ξ(ɛ^) for ɛR ranging between ɛ and ɛions; that is, there is a one-to-one relationship between R and intermediate dielectric constants ɛR. (Note that we have only used the existence of ξ to establish the existence of the ɛR; the exact form of ξ is not required for this derivation nor for the 0 results below that approximate how hD ij (x,x ;R) depends on R.) With this one-to-one relationship in mind, the dielectic TCF 0 hD ij (x,x ;R) is the ion cloud that accumulates at a wall with dielectric constant ɛR because of the surface charge induced by an ion of species i at x and the ions that screen it. To approximate the effective charge of this screened ion, we consider two limits. First, when x is far from the dielectric interface, the screening shell is complete, and the effective charge is 0. Second, when the ion is touching the dielectric interface, the ion’s charge is so close that it dominates the induced charge because the finite size of the ions keeps other ions at least a distance of Rij away; the screened ion’s effective charge is zi. The change in effective charge between these two limits occurs as the screening ions fill in the remaining void between the dielectric interface and the central ion, namely, the cylinder of radius Ri and height |sx x| (sx is the point on S closest to x). Here, we approximate the distance dependence as changing exponentially with the Debye length λ, namely, exp((Ri |sx x|)/λ). While an approximation (which is validated a posteriori in the figures), it is the kind of distance dependence that one would expect by integrating the counter charge from linearized PoissonBoltzmann (LPB) theory in that cylinder. The effective surface charge then is ER Þ σeff i ðs; x, ^
ðRi jsx xjÞ=λ
¼ σi ðs; x, ^ER Þe
ð9Þ
because σi is proportional to the ion’s effective charge (eq 1). Next, to examine how ions accumulate around such a surface charge, we consider the screened potential around each surface element s. Because each individual surface element is small, we approximate this potential as a Yukawa potential (like in LPB theory) so that γ ðxÞ σeff ER Þ js x0 j=λ i ðs; x, ^ e φi ðs, x ; x, ^ER Þ i 0 4πE0 js x j 0
ð10Þ
with an as-yet-undetermined coefficient γi(x). Because the charge on a single surface element is infinitesimal, the potentials from all of the elements are (approximately) additive. Therefore, the potential from the entire surface charge at x0 is Z Φi ðx 0 ; x, ^ER Þ φi ðs, x 0 ; x, ^ER Þ ds ð11Þ S
0
With another LPB approximation, for |x x | g Rij zj e ext 0 hDij ðx, x 0 ; RÞ Φi ðx0 ; x, ^ER ÞeVj ðx Þ=kT kT
ð12Þ
and 0 when ions overlap. As before, this TCF is multiplied by the Bolzmann factor of the ions’ accessible volume, and it is rescaled so that the ion cloud is charge-neutral, giving Z 4πEE0 kT σi ðs; x, ^E Þ ds S Z γi ðxÞ ¼ σ E Þ js x0 j=λ i ðs; x, ^ bath 2 e χj ðx0 Þ dx 0 e2 zj Fj 0j js x 0 jx x j g Rij j
Note that while the arrangement of ions around a screened surface charge requires the most approximations, we show their validity a posteriori by comparing the resulting DFT with MC simulations (see the figures later). 0 To understand how hD ij (x,x ;R) depends on R, we again use the existence of the separation σi(s;x,ɛ^R) ≈ ξ(ɛ^R) σi*(s;x). There0 ^R) (by eq 9), so that fore, hD ij (x,x ;R) ξ(ɛ hDij ðx, x 0 ; RÞ RhDij ðx, x 0 Þ
ð13Þ
D 0 0 where hD ij (x,x ) = hij (x,x ;1). (Note that the existence of the factoring of σi(s;x,ɛ^R) is used only to estimate the dependence of D 0 0 hD ij (x,x ;R) on R. In the actual calculation of hij (x,x ), one should use the exact σi(s;x,ɛ^) from eq 1.) The DFT is derived by substituting eqs 48, 12, and 13 into eq 3. Integrating over R gives the functional Z Fi ðxÞψDii ðx, xÞ dx 2F D ½fFk ðxÞg ZZ þ ΔFi ðxÞΔFj ðx 0 ÞHij ðx, x 0 ; 3ÞψDij ðx, x 0 Þ dx dx 0 ZZ 0 0 D 0 0 þ ΔFi ðxÞFÆ0æ j ðx ÞHij ðx, x ; 2Þψij ðx, x Þ dx dx ZZ 0 D 0 0 þ ΔFj ðx0 ÞFÆ0æ i ðxÞHij ðx, x ; 2Þψij ðx, x Þ dx dx ZZ Æ0æ 0 0 D 0 0 þ FÆ0æ ð14Þ i ðxÞFj ðx ÞHij ðx, x ; 1Þψij ðx, x Þ dx dx
where, for space considerations, the Einstein summation convention is used, ΔFk(y) = Fk(y) FÆ0æ k (y), and 1 1 D 0 wall 0 h ðx, x 0 Þ Hij ðx, x0 ; nÞ ¼ ðhion ij ðx, x Þ þ hij ðx Þ þ 1Þ þ n n þ 1 ij ð15Þ For the functional derivative 1 Æ0æ δF D =δFi ðxÞ zi eðφÆ0æ CD ðxÞ φC ðxÞÞ 6 1 1 þ ψDii ðx, xÞ þ zi eðφCD ðxÞ φC ðxÞÞ 2Z 3 1 0 D 0 0 cj ðx 0 ; 2Þhion þ ij ðx, x Þψij ðx, x Þ dx 3 Z 1 0 D 0 0 cj ðx 0 ; 2Þhwall þ ij ðx, x Þψij ðx, x Þ dx 3 Z 1 cj ðx 0 ; 3ÞhDij ðx, x 0 ÞψDij ðx, x 0 Þ dx0 þ 4
ð16Þ
where 1 cj ðx 0 ; nÞ ¼ Fj ðx0 Þ þ FÆ0æ ðx 0 Þ n j
ð17Þ
0 This used the fact that, for ψij(x,x0 ) = ψCij (x,x0 ) þ ψD ij (x,x ) Z Fj ðx 0 Þψij ðx, x 0 Þ dx 0 ¼ zi eφCD ðxÞ ð18Þ
∑j
∑
∑j 1180
Z
Fj ðx 0 ÞψCij ðx, x 0 Þ dx0 ¼ zi eφC ðxÞ
ð19Þ
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER 19
where φC and φCD satisfy the Poisson equations E0 r 3 ½EðxÞrφCD ðxÞ ¼ e E0 Er2 φC ðxÞ ¼ e
∑j zj Fj ðxÞ
∑j zj Fj ðxÞ
ð20Þ ð21Þ
Æ0æ Æ0æ Æ0æ Similar definitions hold for φC and φCD using Fj (x0 ). Next, we compare the DFT calculations to MC simulations. In general, the DFT took ∼1 min on a single Macintosh PowerPC G5 2 GHz processor, which is substantially faster than the ∼1 h for the MC but slower than the Coulombic DFT, which uses fast Fourier transforms to compute convolutions created by the purely radial interaction potential.20 Figures 13 show the comparison of the ARMY dielectric functional at a single planar dielectric interface to MC simulations (done as previously described6). To focus solely on the effect of the dielectric interface, there is no additional fixed surface charge on the wall.
Figure 1. Concentration profiles of ions near a single planar dielectric interface at x = 0; ɛ is 2 (black), 20 (red), 40 (green), 80 (blue), and ¥ (magenta). Monovalent ions (diameter 3 Å) at bath concentrations of 0.01, 0.1, and 1 M. The dashed line shows the results of using only the self-energy term for ɛ = 2. In this and all figures, the lines are the DFT results, and the symbols are from MC simulations; the surface charge is 0, and T = 298.15 K, ɛions = 80, and ɛ = 500 were used in the DFT to approximate ¥.
Even this seemingly simple system is challenging because dielectric interactions are dominated by correlations. This can be seen by considering the restricted primitive model (RPM, where all ions have the same size) of monovalent ions. In this case, both cations and anions feel the same potential ψD ij because ψD ij zizj, and therefore, the cation and anion concentration profiles are identical and charge-neutral everywhere (Figure 1). If the correlations among ions are ignored in a mean-field approach, the mean-field averaged potential is identically 0 because the mean electrostatic potential φCD(x) = 0 everywhere, and therefore, δFD/δFi(x) = ψD ii (x,x)/2. This approximation is very poor, especially at high concentrations, as shown for ɛ = 2 with dotted lines in Figure 1. In this geometry, 1 19 and it is the hion ψD ii (x,x) decays as x , ij term in eq 16 that compensates for this far from the wall (not shown). This makes sense intuitively because this correlation describes the ion cloud around the central ion that produced the self-energy term. The hD ij term in eq 16 is important near the interface, sometimes becoming as large as the hion ij term (not shown). The contributerm in eq 16 was relatively small in all of the tion of the hwall ij cases studied (not shown) but will be larger when a fixed surface charge is included. The relative sizes of these terms then give some physical insight into the repulsion of ions by low dielectrics; the asympotic behavior is governed by the canceling of the self-energy and bulk screening cloud terms, while near the dielectric interface (∼8 Å for Figure 2A), the ion cloud near the induced surface charge and the distorted screening cloud around each ion dominate the formation of the double layer. The TCFs described above are based on principles derived from LPB, including additivity and eq 12. Therefore, it must be tested whether the functional fails when these approximations are questionable, namely, high valence and high concentration. This is done in Figures 2 and 3 for ion sizes and concentrations of normal electrolytes. Further testing of other conditions (e.g., significantly higher concentrations of ionic liquids5) will be reported in future work. Figure 2 shows the RPM for divalent cations and monovalent anions. Over the range of ɛ = 2¥ for = 0.5Fbath = 1 M, the ARMY dielectric functional works Fbath 2 1 extremely well. In fact, it reproduces many subtle features of the concentration profile found in the MC simulations. These include the steep initial slope (near x j 3 Å) in the anion profile, which then decreases, something not seen in the cation profile.
Figure 2. (A) Concentration profiles of divalent cations and monovalent anions (diameter 3 Å) at a 1 M bath cation concentration for ɛ = 2 (black), 20 (red), 40 (green), 80 (blue), and ¥ (magenta). (B) ɛ = ¥ for the same ions at 0.01 (black), 0.1 (blue), and 1 M (magenta) bath cation concentrations. Solid lines/shaded symbols: cations. Dashed lines/open symbols: anions. 1181
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
The Journal of Physical Chemistry Letters
LETTER
’ ACKNOWLEDGMENT This material is based on work supported by, or in part by, the U.S. Army Research Laboratory and the U.S. Army Research Office under Contract No. W911NF-09-1-0488. Special thanks to Dezs€o Boda for providing the Monte Carlo simulations. Very useful discussions with Roland Roth and Jaydeep Bardhan are gratefully acknowledged. ’ REFERENCES
Figure 3. Concentration profiles of Kþ (solid lines/shaded symbols, diameter 2.76 Å), Ca2þ (dashed lines/open symbols, diameter 2 Å), and Cl (dotted lines/half-filled symbols, diameter 3.62 Å) at bath concentrations of 1, 0.1, and 1.2 M, respectively. The main part shows ɛ = ¥, and the inset shows ɛ = 2.
This is another physical insight: the monovalent ions actually have a shaper depletion near the dielectric interface than the divalents. The DFT also reproduces the small “bump” in the cation profile, except that the maximum of this bump is shifted to larger x. This and the tendencies for the DFT to be above the MC near the interface and below thereafter are not, however, defects of the dielectric functional. They mirror the same features in the ion profiles when there is no dielectric difference between the bath and the wall (ɛ = 80, blue lines). This is true for the vast majority of the profiles studied so far, indicating that many of the errors seen in the curves with dielectrics are due to the underlying Coulomb electrostatics functional. Figure 3 shows the results for ions of different size, specifically a mixture of Kþ, Ca2þ, and Cl . The main part of Figure 3 shows the results for ɛ = ¥. While the DFT tends to overestimate the adsorption of all ion species (a pattern also seen in Figure 1), the agreement with MC is generally good, with DFT reproducing the minimum in the Kþ profile and Ca2þ being adsorbed more than Kþ. The inset shows very good agreement for ɛ = 2. The “kink” in the Cl curve for ɛ = ¥ is an artifact for the largest ion only, stemming from the multiplication of bulk TCF by the Boltzmann factor of the wall geometry (data not shown). Its effect is small in the cases considered so far, but it is an area of improvement for the ARMY approach. Dielectric interfaces have first-order effects on the structure of ions that cannot be reproduced by mean-field approaches. The DFT derived here includes the necessary correlations to better understand the dielectric electrical double layer over a wide range of dielectric constants, ion concentrations, and ion species (including mixtures). In this first report, the planar geometry is considered, but the theory is general and will be tested in other geometries in future work. Even in the planar geometry, it is the first theory of ions with different sizes at a dielectric interface. Moreover, this work shows that a dielectric theory can be constructed as a perturbation of the Coulomb potential.
’ AUTHOR INFORMATION Corresponding Author
*E-mail: [email protected].
(1) Boda, D.; Fawcett, W. R.; Henderson, D.; Sokolowski, S. Monte Carlo, Density Functional Theory, and PoissonBoltzmann Theory Study of the Structure of an Electrolyte Near an Electrode. J. Chem. Phys. 2002, 116, 7170–7176. (2) \Boda, D.; Valisko, M.; Eisenberg, B.; Nonner, W.; Henderson, D.; Gillespie, D. Combined Effect of Pore Radius and Protein Dielectric Coefficient on the Selectivity of a Calcium Channel. Phys. Rev. Lett. 2007, 98, 168102. (3) van der Heyden, F. H. J.; Stein, D.; Besteman, K.; Lemay, S. G.; Dekker, C. Charge Inversion at High Ionic Strength Studied by Streaming Currents. Phys. Rev. Lett. 2006, 96, 224502. (4) Torrie, G. M.; Valleau, J. P. Double Layer Structure at the Interface Between Two Immiscible Electrolyte Solutions. J. Electroanal. Chem. 1986, 206, 69–79. (5) Loth, M. S.; Skinner, B.; Shklovskii, B. I. Anomalously Large Capacitance of an Ionic Liquid Described by the Restricted Primitive Model. Phys. Rev. E 2010, 82, 056102. (6) Boda, D.; Gillespie, D.; Nonner, W.; Henderson, D.; Eisenberg, B. Computing Induced Charges in Inhomogeneous Dielectric Media: Application in a Monte Carlo Simulation of Complex Ionic Systems. Phys. Rev. E 2004, 69, 046702. (7) Buyukdagli, S.; Manghi, M.; Palmeri, J. Variational Approach for Electrolyte Solutions: From Dielectric Interfaces to Charged Nanopores. Phys. Rev. E 2010, 81, 041601. (8) Hatlo, M. M.; Lue, L. Electrostatic Interactions of Charged Bodies from the Weak- to the Strong-Coupling Regime. Europhys. Lett. 2010, 89, 25002. (9) Outhwaite, C. W.; Bhuiyan, L. B. An Improved Modified PoissonBoltzmann Equation in Electric-Double-Layer Theory. J. Chem. Soc., Faraday Trans. 2 1983, 79, 707–718. (10) Evans, R. Density Functionals in the Theory of Nonuniform Fluids. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; pp 85176. (11) Tang, Z.; Scriven, L. E.; Davis, H. T. Interactions Between Primitive Electrical Double Layers. J. Chem. Phys. 1992, 97, 9258–9266. (12) Kierlik, E.; Rosinberg, M. L. Density-Functional Theory for Inhomogeneous Fluids: Adsorption of Binary Mixtures. Phys. Rev. A 1991, 44, 5025–5037. (13) Rosenfeld, Y. Free Enery Model for Inhomogeneous Fluid Mixtures: Yukawa-Charged Hard Spheres, General Interactions, and Plasmas. J. Chem. Phys. 1993, 98, 8126–8148. (14) Gillespie, D.; Nonner, W.; Eisenberg, R. S. Coupling Poisson NernstPlanck and Density Functional Theory to Calculate Ion Flux. J. Phys.: Condens. Matter 2002, 14, 12129. (15) Gillespie, D.; Valisko, M.; Boda, D. Density Functional Theory of the Electrical Double Layer: The RFD Functional. J. Phys.: Condens. Matter 2005, 17, 6609–6626. (16) Bardhan, J. P.; Knepley, M. G.; Anitescu, M. Bounding the Electrostatic Free Energies Associated with Linear Continuum Models of Molecular Solvation. J. Chem. Phys. 2009, 130, 104108. (17) Sullivan, D. E. Statistical Mechanics of a Nonuniform Fluid with Long-Range Attractions. Phys. Rev. A 1982, 25, 1669. (18) Davis, H. T. Statistical Mechanics of Phases, Interfaces, and Thin Films; Wiley-VCH: New York, 1996. (19) Jackson, J. D. Classical Electrodynamics, 3rd ed.; John Wiley and Sons, Inc.: New York, 1999. (20) Knepley, M. G.; Karpeev, D.; Davidovits, S.; Eisenberg, R. S.; Gillespie, D. An Efficient Algorithm for Classical Density Functional Theory in Three Dimensions: Ionic Solutions. J. Chem. Phys. 2010, 132, 124101. 1182
dx.doi.org/10.1021/jz2001908 |J. Phys. Chem. Lett. 2011, 2, 1178–1182
