Atomistic and Molecular Effects in Electric Double Layers at High

Jun 16, 2015 - Mechanical Engineering Department, Stanford University, Stanford, ... Multiscale Model for Electrokinetic Transport in Networks of Pore...
0 downloads 3 Views 2MB Size
Download PDF
Article pubs.acs.org/Langmuir

Atomistic and Molecular Effects in Electric Double Layers at High Surface Charges Jonathan W. Lee,† Ali Mani,‡ and Jeremy A. Templeton*,† †

Thermal/Fluid Science & Engineering Department, Sandia National Laboratories, Livermore, California 94550, United States Mechanical Engineering Department, Stanford University, Stanford, California 94305, United States

‡

ABSTRACT: The Poisson−Boltzmann theory for electrolytes near a charged surface is known to be invalid due to unaccounted physics associated with high ion concentration regimes. To investigate this regime, fluids density functional theory (f-DFT) and molecular dynamics (MD) simulations were used to determine electric surface potential as a function of surface charge. Based on these detailed computations, for electrolytes with nonpolar solvent, the surface potential is shown to depend quadratically on the surface charge in the high charge limit. We demonstrate that modified Poisson−Boltzmann theories can model this limit if they are augmented with atomic packing densities provided by MD. However, when the solvent is a highly polar molecule, water in this case, an intermediate regime is identified in which a constant capacitance is realized. Simulation results demonstrate the mechanism underlying this regime, and for the salt water system studied here, it persists throughout the range of physically realistic surface charge densities so the potential’s quadratic surface charge dependence is not obtained.

■

theory which prevents ion concentrations from exceeding 1/d3, where d is the ion diameter. In the high charge limit, the model predicts an electric surface potential of

INTRODUCTION With the increased interest in electrical energy storage, research on the interface between electrolytes and electrified substrates at high voltages has increased. Various theories for the electric double layer (EDL) have been developed and studied, most notably mean-field theories such as Poisson−Boltzmann/ Gouy−Chapman theory and subsequent modifications to include steric effects.1−10 The original Gouy−Chapman model assumes a Boltzmann distribution of ions sufficient to screen the surface charge with a characteristic length scale referred to as the Debye length, λD = [(ϵkBT)/(2noe2)]1/2, where ϵ is the permittivity, kB is the Boltzmann constant, T is the temperature, no is the bulk molarity, and e is the elementary charge (multiplied by the ion valency). This model allows for unlimited ion concentrations near the solid−fluid interface, surpassing physically allowable close-packing density for surface charge magnitudes common to microfluidic pumps, electrochemical sensors, and supercapacitors. In general, the permittivity varies depending on the local medium, i.e., ϵ = ϵr ϵ0, where ϵr is the medium’s relative permittivity, and ϵ0 is the vacuum permittivity. The main assumptions of most Gouy−Chapman-type models are (i) ϵ is constant, and (ii) the only significant interactions are Coulombic.1−3 A third assumption (which is implied by the two but worth noting) is that the solvent does not affect the equilibrium configuration of the system, aside from providing a dielectric medium. Typical modifications to the Gouy−Chapman theory involve statistical mechanical derivations to include a Stern layer3 immediately adjacent to the surface. Borukhov et al.,7 for example, presented a modified Poisson−Boltzmann (MPB) © XXXX American Chemical Society

ζe 1 d3σ 2 = kBT 2 ϵkBT

(1)

where σ is the surface charge density. Typical surface charges may vary depending on electrode material, but they are normally bounded by 6(1) C/m2. Several of the assumptions regarding the dielectric medium have also been relaxed in PB theory through the addition of Langevin dipoles by Iglic̆ et al.11 These dipoles act to both provide a variable relative dielectric coefficient more representative of a polar solvent as well as displace ions near the surface. Kilic et al.8 developed a composite diffuse layer (CDL) model assuming a “condensed layer” in place of the Stern layer. The condensed layer is a continuum region of counterions with uniform density (co-ions are excluded from the condensed layer). At the interface between the condensed layer and the outer diffuse layer, the electric field is used for boundary condition matching. The CDL model simplified for the high charge limit yields the same behavior as the MPB model. A similar approach to account for atomic packing was introduced by Bohinc et al.12 In addition to Coulombic interaction between ions, they added interactions via a Yukawa potential to the free energy and derived a Poisson-Helmholtz-Bolzmann Received: January 23, 2015 Revised: April 27, 2015

A

DOI: 10.1021/acs.langmuir.5b00215 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

molecular counterions. In the first system, we model two weakly interacting (ideally noninteracting), oppositely charged EDLs and utilize the ionic antisymmetry to simultaneously collect data for different ion sizes. In order to examine more realistic systems, the second system features a monovalent tetrahedral ion (e.g., tetrafluoroborate) as the counterion with a graphite electrode. The co-ion is the same generic spherical ion as in the previous system, and here we only model a single EDL with an approximate bulk electrolyte in the far-field. Additional information can be found on these modeling strategies in Lee et al.16 and the references therein. Here, we employ the Tramonto software package17 to perform our f-DFT calculations. Details on the numerical methods implemented in this software can be found in Heroux et al.18 Atoms interact electrostatically as well as via LJ. Given an LJ diameter d, the inner cutoff dhs is 0.93 d (this fraction is a common choice because it gives the best comparison between f-DFT, MC, and MD simulations16), and the outer cutoff is 4.127 d. Inside of 0.93 d, the fluid atom is treated as a hard sphere. The hard sphere functional is Rosenfeld’s fundamental measure theory.19 Electrostatics are handled through solving Poisson’s Equation with the mean spherical approximation, as described in Lee et al.16 The f-DFT domain is 15 atomic diameters long, where the solid occupies the first full atomic diameter (i.e., the center of the first fluid atom is located at approximately 1.5 atomic diameters from the domain edge). A structureless substrate provides one boundary condition and is modeled with an LJ 10-4-3 potential20 of the form

equation governing the EDL. Their work also generalizes to other types of ion−ion interactions. To the best of our knowledge, these models have not been quantitatively verified at the high charge limit against experiments or a high fidelity numerical simulation. In this work, we use fluids density functional theory (f-DFT) and molecular dynamics (MD) simulations to rigorously investigate the charge-potential relations without resorting to a macroscopic model. Specifically, we explore the upper limit of surface charge density and its effect on ionic packing and capacitance. The density functional theoretic approach has been demonstrated to account for ion size effects in the EDL by Gillespie et al.13 and be one of the more explanatory theories for EDL capacitance and ion distributions14 as compared to Monte Carlo simulations. Therefore, we will explore an algebraic simplification of MPB theory, the higher-fidelity but computationally efficient f-DFT, and the more expensive but more accurate MD methods for the EDL. Our results show the details of the precise close-packing arrangement in such a limit, which confirms these theories if the atomic packing structure identified by modeling is included. We further demonstrate that the presence of solvent in these layers significantly alters the predicted behavior due to changes in the electrical medium and nonconstant effective volumes of the solvent.

■

NONPOLAR SOLVENT RESULTS We first compare simulations of atomic solutions with a background uniform permittivity as assumed by most existing theories. For f-DFT, we simulate a three-component model (3CM) fluid near a structureless solid substrate. The three atomic components of the 3CM are two ionic species (positively and negatively charged) and a charge-neutral solvent atom.15 Counterions with varying valency are tested. The permittivity of the medium is handled through an effective uniform dielectric constant. In our 3CM, ions and solvent atoms all have the same Lennard-Jones (LJ) parameters (see Table 1). Due to the symmetry, f-DFT calculations need only

⎤ ⎡ 2 d10 d4 2 d3 v(z) = 2π ϵLJ⎢ 10 − 4 − 3⎥ ⎣5 z z 3(z + (0.61/ 2 )d) ⎦ (2)

where z is the wall-normal position measured from the wall. Outside of 3.5d, the potential is cut off and shifted by subtracting the value at the cutoff. A relative permittivity value of ϵr = 80 is applied to mimic the polarization of water. Simulations are performed at room temperature of 300 K. The external potential is equivalent to an electric field of magnitude σ/ϵ. These parameters are all listed in Table 1. Atomic number density profiles can be readily extracted from f-DFT and MD simulations. An example of this can be seen in Figure 1a. Scaling atomic densities by their valence yields the free charge density (Figure 1b). Gauss’ Law and Poisson’s Equation can be used to determine electric field (Figure 1c) and electric potential (Figure 1d) profiles, respectively. Oscillations arising in these plots are due to layers of atoms undergoing thermal motions separated by space devoid of particles. MPB theories are predicated on a mean-field representation which typically neglects the specific atomic details and effectively averages over the layers. Atomic-level calculations cannot therefore be directly spatially compared to mean-field theories, but integrated quantities are expected to be consistent between both representations. For the purpose of comparing with theory, the “electric surface potential” is defined to be the electric potential measured at the fluid−solid interface, which is required to compute and compare integrated quantities such as the electric potential drops. The fluid−solid interface is taken to be one-half atomic diameter away from the LJ 10-4-3 potential’s origin, consistent with Lee et al.16 In this work, we examine the total potential drop across the EDL rather than the zeta potential, because in the high surface

Table 1. Simulation Parameters for f-DFT and MD Calculationsa parameter

f-DFT

MD (1)

Atom−Atom ϵ (kcal/mol) Atom−Atom d (Å) Atom−Atom inner cutoff (×d) Atom−Atom outer cutoff (×d) Atom−Wall 2πϵLJ (kcal/mol) Atom−Wall d (×d) Atom−Wall outer cutoff (×d) Relative Permittivity ϵr Nominal Temperature (K) Atomic Mass (g/mol)

0.4958 3.1507 0.93 4.127 0.4958 1.0 11.0275 80 300 N/A

0.4958 3.1507 N/A 4.127 0.4958 1.0 11.0275 80 300 18.0154

LJ

MD (2)

N/A 4.127 0.4958 1.0 11.0275 80 300

a

MD (1) refers to the generic spherical ions analogous to the f-DFT calculations. MD (2) refers to the cases with tetrahedral monovalent counterions. The word “atom” is used to represent either a solvent or solute particle, as opposed to the LJ wall.

model one EDL, with the far-field bounded by bulk conditions. Bulk conditions are counterion concentrations of 1 M and a bulk solvent density, ρ* = ρ d3, of 0.7 (where d is the atomic diameter). For MD, we simulate an analogous system of generic spherical (monovalent only) ions in addition to a system with B

DOI: 10.1021/acs.langmuir.5b00215 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Figure 2. Charge-potential relation comparison of f-DFT and MD results. The theoretical model prediction (dashed line) compares remarkably well at high surface charges.

estimate of the upper limit to the surface charge density. However, the molecular counterion simulations were performed up to a (dimensional) surface charge of approximately 0.73 C/m2. To nondimensionalize the latter simulation’s results, we approximate the tetrahedron as a sphere, such that the radius is equivalent to the bond length plus half the vertex atom’s LJ diameter (i.e., the distance from centroid to “edge” of vertex atom). In general, for larger charged molecules, e.g., ionic liquids, this theory will become an increasingly useful approximation. As shown in Figure 2, the f-DFT results (circles) and MD results (diamonds) approximately follow the same trend. Minimum and maximum potentials observed in the course of the MD simulations are denoted by the error bars to estimate the uncertainty in those measurements. The parabolic shape is indicative of a potential created by a uniform density of charges. Here, we present a simplified model by focusing only on the packing effects and neglecting all other physics associated with steric effects and the diffuse layer. As shown in Figure 3, we Figure 1. Example of f-DFT/MD EDL comparison. Electrolyte concentration is targeted to 1 M, and bulk solvent density, normalized by d3, is targeted to 0.7. This example corresonds to a σ* = 3.0. MD electrolyte concentration is 1.02 M, and bulk solvent density is 0.704.

charge limit, only a vanishingly small fraction of the potential drop occurs in the diffuse region. As predicted by MPB theories and borne out in the simulations, most of the potential drop takes place over the condensed layers of atoms. Figure 2 shows the electric surface potential as a function of applied surface charge. The inverse of its slope is the differential capacitance of the EDL. The left-hand figure is plotted on a linear scale, whereas the right-hand figure is on a log−log scale. The f-DFT calculations were performed with normalized surface charges (σ* = σ d2/e) up to 6.5. Convergence for higher charge cases was difficult to achieve, given the strong ordering within the compact layer. Depending on the ion and solvent sizes, eventually the dimensionless surface charge density σ* will become unphysical as chemical processes not included in the MD and f-DFT models become controlling. For an effective salt-water model where the distance between constituents is roughly the same size as the spacing of atoms in a metal lattice, σ* ≈ 1 corresponds to a dimensional surface charge density of approximately 1 C/m2. This value is a good

Figure 3. (a) Schematic of charge distribution in an electrolyte next to a highly charged surface as proposed by Kilic et al.,8 and (b) our idealized model consisting of only counterions with maximum packing and total charge equal to the surface charge.

assume that in the limit of large surface charge, the potential drop in the EDL is governed by a region of densely packed multilayered counterions, i.e., a condensed layer. Figure 4 shows example snapshots from our MD simulations. The crosssectional view (Figure 4b) reveals a triangular packing arrangement, which we infer to be to a hexagonal close-packed or face-centered cubic lattice. Both lattices yield the maximum C

DOI: 10.1021/acs.langmuir.5b00215 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

integrating eq 5. Given the conditions in Table 1, eq 6 reduces to 2 ζe e 2 ⎛ σd 2 ⎞ = 0.28 ⎜ ⎟ kBT d ϵkBT ⎝ e ⎠

(7)

By putting eq 1 in a similar form, the difference between the standard high charge limit theories and the present simulations can be understood 2 1 e 2 ⎛ σd 2 ⎞ ζe = ⎟ ⎜ kBT 2 d ϵkBT ⎝ e ⎠

The difference in the two equations is the prefactor which accounts for the volume occupied by ions. The prefactor in eq 7 is approximately 60% of the MPB7 and CDL8 predictions. Equation 7 is plotted as a dashed line in Figure 2. The model appears to match simulation data remarkably well on the linear scale. On the log−log scale, discrepancies are emphasized in the low charge regime due to ion and solvent layering/intercalation observed in the simulations. It is in this low charge regime where PB theory is observed to be congruent with simulation results.16 In contrast to studies at lower charges, in Figure 2, as well as subsequent figures except when noted, the total potential drop across the double layer is plotted rather than the zeta potential. For many applications, particularly electrokinetic flows, the zeta potential is more relevant because it is related to the number of mobile ions as well as the electrical force they experience. In the present high-charge regime, nearly all of the screening ions are immobile and the potential drop outside of the condensed layer is very small, as shown, e.g., in Figure 1. This regime is more relevant to electrical energy storage devices. It should be noted that strong oscillations near the walls are observed in the density profiles from both types of simulation. These lead to step-like functions in the electric field plots and inflection points in the electric potential plots. Those inflection points correspond to inflection points in the charging curve (e.g., Figure 7 in Lee et al.16). As shown in Lee et al., the capacitance in between each of those inflection points is approximately constant with new layers becoming active for each unit increment of dimensionless surface charge. Conversely, the presented continuum model assumes that the counterion density in the near-wall region is a uniformly distributed block of charge. Despite this simplification, the trends are still observed to follow the same behavior and be quantitatively well-predicted because of capturing the most important atomistic effect, namely, packing density, and incorporating it into the model. This indicates that in the high charge regime, density oscillations arising from atomic packing are not particularly important. In other words, for the high surface charge regime, the inflection points in the capacitance trend appear so closely spaced that the curve is essentially quadratic. The diffuse layer’s contribution to screening is negligible in this limit.

Figure 4. Snapshots from our MD simulations. (a) Side view of ionic packing arrangement. (b) Cross-sectional view of one ionic layer near the surface.

packing fraction of 0.74. The thickness of this condensed layer is assumed to be such that the total charge balances the applied surface charge. In this model, we simply ignore contributions of the diffuse layer to the potential drop. As we shall see, this approximation is valid since the diffuse layer can at most contribute potential drops on the order of a few thermal voltages, which is negligible compared to the potential across the condensed layer. As such, the homogenized charge density of the compact layer can be estimated as ρmax =

0.74e 6 × 0.74 e = 3 π 4/3π(dhs/2)3 dhs

(3)

Assuming a surface charge density of σ, the thickness of the condensed layer that would balance the wall charge is σ λc = ρmax (4) This double layer thickness is only a valid measure in the high charge limit when nearly the entirety of the potential drop occurs across the condensed layer. In that sense, it is the counterpart to the Debye length arising from PB theory in the low charge limit. As both thickness estimates arise in limiting conditions, more general length scales are required at intermediate surface charge regimes, such as the screening fraction definition of Bohinc et al.,21 but the present work focuses solely on the high-charge limit. Variation of electric potential, ϕ, across the double layer can be described using Poisson’s Equation

ρ d2ϕ = − max 2 ϵ dz

(5)

Upon integrating twice and matching the field on the surface (dϕ/dz = σ/e) we obtain a relation between the surface charge and the electric surface potential, ζ, i.e., the potential drop across the double layer. Using eqs 3 and 4, the result can be written in the following dimensionless form 2 ⎡ ⎛ dhs ⎞3⎤ e 2 ⎛ σd 2 ⎞ ζe = ⎢0.3538⎜ ⎟ ⎥ ⎟ ⎜ ⎝ d ⎠ ⎥⎦ d ϵkBT ⎝ e ⎠ kBT ⎣⎢

(8)

■

POLAR SOLVENT RESULTS Having established the viability of existing MPB theories for simple nonelectrostatically active solvents, we turn our attention to the effect of polarized solvent molecules. In this case, to eliminate the effects of differing ion sizes, a uniform ion is considered in conjunction with the TIP3 model of water22,23 such that the only asymmetry in the system will be due to

(6)

In eq 6, the term in the brackets is the dimensionless maximum packing density, while the term multiplying it arises from twice D

DOI: 10.1021/acs.langmuir.5b00215 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

charges considered in this part of the study. This is in contrast to the results of MPB theories which predicted a quadratic behavior in the high charge limit. Note that because the reference for the zero location varies slightly when water is the solvent, the comparison in Figure 5 is only approximate. Also of note is that the potential drops are different at the anode and the cathode, while MPB theories predict no difference. The difference is due to different packing structures arising from hydrogen-ion alignment for negative ions and hydrogen-ion offset for positive ions. As a result, charges can be packed more efficiently at the anode. The presence of this linear regime can be understood by considering the behavior of filling the adsorbed layers. However, what is different in this case is the surface charge ranges over which the filling process persists. MPB theories predicted a faster layer-by-layer filling and thus linear growth of double-layer thickness with charge, resulting in a quadratic potential-charge relation. In contrast, the presented molecular solvent models predict formation of a single dominant layer over these surface charge ranges (MPB theories predict 2−3) which continuously fills, and as such results in an approximately linear potential-charge relation. Of particular note is that MPB theories use the bulk permittivity adjacent to the wall, which is not the case with the explicit solvent model, so the persistence of the observed deviation over such a range of surface charges is even more significant. To investigate the causes of this phenomena, layer-specific quantities are computed by integrating over a layer of particles, which is identified by minima in density profiles. We restrict our attention to the first layer as it reveals the reason behind this behavior and layer integrals can be more precisely computed. Figure 6 shows the number of counterions present

orientation of the water molecules. As capturing this type of molecular structure in solution is challenging for density functional theory, only MD models are considered following those studied in Lee et al.24 Fourteen simulation sets of five replicas each were performed varying the dimensionless surface charge from 0 to 2.4. For all cases, a target density of 0.95 g/ cm3 and target molarity of 1 mol/L were chosen, with actual densities ranging from 0.94 to 0.96 g/cm3 and actual molarities from 0.84 to 1.17 mol/L. Although we note that for this particular system values beyond 1 mol/L would be difficult to realistically obtain, they are included so that we can fully illustrate an intermediate behavior between the low and high charge limits which we anticipate to be increasingly relevant to systems with larger molecules. Before proceeding with the analysis, it is worth assessing what the expected results would be according to the MPB theories. If the only effect of the solvent were that of a bulk dielectric, then the capacitance of this double layer would be slightly lower than in the nonpolar solvent case, as the relative permittivity for this water model is approximately 72.25 Over this range of surface charges, there would be two inflection points corresponding to the first two adsorbed layers forming. This would be sufficient to observe the quadratic trend in the MPB model. The MD model for this system can be found in Lee et al.24 The important point to recall is that the model for water is nonpolarizable; hence, all changes in polarization are due to rotation of the molecule. To analyze the results, we use postprocessing coarse-graining methods26 for atomic densities as well as molecular coarse-graining theory25 to compute polarization densities throughout the system. Representative density, charge density, electric potential, and polarization profiles can be found in Lee et al.24 and Mandadapu et al.25 Defining the “solid-fluid” interface layer is difficult in this case because the particles increasingly approach the surface with increasing charge density. In contrast to the nonpolar solvent model, we take the reference plane for the potential drop to be at the beginning of the adsorbed layer, i.e., where the density of ions or solvent first becomes nonzero. This quantity can be easily identified as the point when the electric field becomes nonconstant. Figure 5 shows that the potential drop versus surface charge has a linear relationship at the higher surface

Figure 6. Counterion density in the first adsorbed layer in the salt water MD simulations.

in the first layer using the number of ions per unit area made dimensionless using the LJ length scale. As σ* increases, the plot shows an asymmetric behavior between the cathode and the anode, but at higher surface charge a symmetric linear trend is observed. The low charge asymmetry is likely due to the increased polarization on the cathode side, shown in Figure 7, such that the water molecules are more energetically favorable to be in the vicinity of the charged surface. The definition of polarization is taken from eq 15 in Mandadapu et al.25 to be

Figure 5. Charge-potential relation for the salt water obtained from MD simulations. E

DOI: 10.1021/acs.langmuir.5b00215 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

count can be verified to be nonconstant. Only by the highest surface charge in this study is this trend possibly changing. The reason for this can be seen in Figure 9 in which the polarization

Figure 7. Total polarization in the first adsorbed layer in the salt water MD simulations. N

p=

∑ qi(x i − xc) i=1

Figure 9. Total polarization per water molecule present in the first adsorbed layer in the salt water MD simulations.

(9)

where qi is the charge associated with the ith atom of a molecule, xi is its position, and xc is the molecule’s center of mass. In all simulations beyond the establishment of this trend, the first layer accounts for the significant majority of the screening charge, and hence it is this linear filling of the first layer which strongly contributes to the linear relationship. Importantly, even by the highest surface charge in this study, the first layer is far from completely filled by counterions; see the water molecule plot in Figure 8.

per water molecule is observed to asymptote to its minimum value (for this model, −0.11) corresponding to the maximum alignment of the water molecules with an electric field vector. It is not until the water molecules are fully aligned with the surface normal that replacement of water by counterions is observed to occur. This behavior will eventually lead to the high-charge limit MPB theoretical results, albeit at a substantially slower rate due to the energetic favorability provided by the water polarization (with slight adjustments for different ion and molecule effective sizes). The data therefore demonstrate that an intermediate regime exists between the established low (Debye−Hückel) and high (MPB) limits in which adsorbed solvent molecules align with the electric field. Important physics not captured in MPB theories include a higher solvent concentration and changes in the adsorbed lattice structure due to changes in the in-plane solvent molecule size arising from these orientations. In this regime a constant capacitance is observed due to linear variations in counterion concentration with surface charge occurring by (1) reduction in effective molecular size and (2) displacement of solvent molecules by counterions. We note that, for water, it will remain in this regime without approaching the high-charge limit for physically realizable systems, although for solutions with larger ions and/or solvent molecules this need not be the case.

■

Figure 8. Water molecule density in the first adsorbed layer in the salt water MD simulations.

CONCLUSIONS In summary, we used f-DFT and MD simulations to numerically determine the charge-potential relationship in the intermediate- and high-charge regimes for nonpolar solvents. Given a close-packed arrangement of ions near the surface, we confirmed the validity of existing MPB theories if the condensed layer packing structure is correctly accounted for. We model the charge density near the wall as a uniform block of charge sufficient to entirely screen the applied surface charge. The simulations and theoretical model both predict a quadratic dependence of potential on surface charge density, consistent with earlier models. Moreover, the maximum packing density

Water molecules, being polarizable, can have a favorable energy of attraction to the charged surface, which is not captured in the MPB theories. As such, their presence continues at higher surface charges than would otherwise be predicted. Also playing a role is that the effective size of the water molecules in the lattice-like adsorbed layers changes with orientation. This is clearly seen in Figure 8 because the number of water molecules remains relatively constant over charge densities for which the counterion count is linearly increasing and co-ions are completely excluded, and the total particle F

DOI: 10.1021/acs.langmuir.5b00215 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

(9) Fedorov, M.; Kornyshev, A. Towards understanding the structure and capacitance of electrical double layer in ionic liquids. Electrochim. Acta 2008, 53, 6835−6840. (10) Wang, H.; Pilon, L. Accurate simulations of electric double layer capacitance of ultramicroelectrodes. J. Phys. Chem. C 2011, 115, 16711−16719. (11) Iglic̆, A.; Gongadze, E.; Bohinc, K. Excluded volume effect and orientation ordering near charged surface in solution of ions and Langevin dipoles. Bioelectrochemistry 2010, 79, 223−227. (12) Bohinc, K.; Shrestha, A.; May, S. The Poisson-HelmholtzBoltzmann model. Eur. Phys. J. E 2011, 34, 1−10. (13) Gillespie, D.; Khair, A.; Bardhan, J.; Pennathur, S. Efficiently accounting for ion correlations in electrokinetic nanofluidic devices using density functional theory. J. Colloid Interface Sci. 2011, 359, 520− 529. (14) Gillespie, D. A review of steric interactions of ions: Why some theories succeed and others fail to account for ion size. Microfluidics and Nanofluidics 2014, 18 (5−6), 717−738. (15) Tang, Z.; Scriven, L.; Davis, H. A three-component model of the electrical double layer. J. Chem. Phys. 1992, 97, 494−503. (16) Lee, J.; Nilson, R.; Templeton, J.; Griffiths, S.; Kung, A.; Wong, B. Comparison of Molecular Dynamics with Classical Density Functional and Poisson-Boltzmann Theories of the Electric Double Layer in Nanochannels. J. Chem. Theory Comput. 2012, 8, 2012−2022. (17) Tramonto Software. https://software.sandia.gov/DFTfluids/. (18) Heroux, M.; Salinger, A.; Frink, L. Parallel segregated Schur complement methods for fluid density functional theories. Siam Journal on Scientific Computing 2007, 29, 2059−2077. (19) Rosenfeld, Y. Free-energy model for the inhomogeneous hardsphere fluid mixture and density-functional theory of freezing. Phys. Rev. Lett. 1989, 63, 980−983. (20) Magda, J.; Tirrell, M.; Davis, H. Molecular dynamics of narrow, liquid-filled pores. J. Chem. Phys. 1985, 83, 1888−1901. (21) Bohinc, K.; Kralj-Iglic̆, V.; Iglic̆, A. Thickness of electrical double layer. Effect of ion size. Electrochim. Acta 2001, 46, 3033−3040. (22) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (23) Price, D.; Brooks, C. A modified TIP3P water potential for simulation with Ewald summation. J. Chem. Phys. 2004, 121, 10096− 10103. (24) Lee, J.; Templeton, J.; Mandadapu, K.; Zimmerman, J. Polar solvent effects on electric double layers in nanochannels. J. Chem. Theory Comput. 2013, 9, 3051−3061. (25) Mandadapu, K.; Templeton, J.; Lee, J. Polarization as a field variable from molecular dynamics simulations. J. Chem. Phys. 2013, 139, 054115. (26) Zimmerman, J.; Webb, E., III; Hoyt, J.; Jones, R.; Klein, P.; Bammann, D. Calculation of stress in atomistic simulation. Modell. Simul. Mater. Sci. Eng. 2004, 12, S319−S332.

and hard-sphere diameter played significant roles in achieving the remarkable agreement between simulations and model. When considering a polar solvent approximating water, different phenomena contribute to an extended intermediate regime distinct from the high- and low-charge limits. In this case, solvent molecules are displaced from the adsorbed layer much more slowly than in MPB models because of the energetic favorability of the dipole in an electric field. The water molecules are observed to first change configuration to accommodate more counterions before they achieve their maximum polarization, at which point displacement occurs. As a result, the linear capacitance regime which exists during layer filling in the MPB theory persists over a broad range of surface charge densities. Physical limits of surface charge density suggest that salt water will not achieve the high-charge limit but can achieve the intermediate regime. Solutions comprising larger solvent and ion molecules may eventually realize the high-charge limit, but likely will have an extended intermediate range due to the same physical processes. Improved constitutive models of the effective molecule size and polarizability will be required to predict the behavior of these electrolytes.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors would like to thank Amalie Frischknecht and Kranthi Mandadapu for helpful comments on a draft of this manuscript. Sandia National Laboratories is a multiprogram laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. Funding for this work was provided by the Laboratory Directed Research and Development (LDRD) and Advanced Scientific Computing (ASC) programs at Sandia National Laboratories, and its support is gratefully acknowledged.

■

REFERENCES

(1) Gouy, L. Sur la constitution de la charge electrique a la surface d’un electrolyte. J. Phys. Theor. Appl. 1910, 9, 457−468. (2) Chapman, D. A contribution to the theory of electrocappilarity. Philosophical Magazine Series 6 1913, 25, 475−481. (3) Stern, O. The Stern theory of electrochemical double layer. Elektrochemistry 1924, 30, 508. (4) Eigen, M.; Wicke, E. The thermodynamics of electrolytes at higher concentration. J. Phys. Chem. 1954, 58, 702−714. (5) Kralj-Iglic, V.; Iglic, A. A simple statistical mechanical approach to the free energy of the electric double layer including the excluded volume effect. J. Phys. II 1996, 6, 477−491. (6) Lamperski, S.; Outhwaite, C.; Bhuiyan, L. A modified PoissonBoltzmann analysis of the solvent primitive model electrical double layer. Mol. Phys. 1996, 87, 1049−1061. (7) Borukhov, I.; Andelman, D.; Orland, H. Steric effects in electrolytes: a modified Poisson-Boltzmann Equation. Phys. Rev. Lett. 1997, 79, 435−438. (8) Kilic, M.; Bazant, M.; Ajdari, A. Steric effects in the dynamics of electrolytes at large applied voltages. I. double-layer charging. Phys. Rev. E 2007, 75, 021502. G

DOI: 10.1021/acs.langmuir.5b00215 Langmuir XXXX, XXX, XXX−XXX