Article Cite This: J. Chem. Eng. Data XXXX, XXX, XXX−XXX
pubs.acs.org/jced
Calculation of Electrical Double Layer Potential Profiles in Nanopores from Grand Canonical Monte Carlo Simulations Evan M. Ney,† Chia-Hung Hou,‡ and Patricia Taboada-Serrano*,§ †
Department of Mechanical Engineering, Rochester Institute of Technology, Rochester, New York 14623, United States Graduate Institute of Environmental Engineering, National Taiwan University, Taipei City, Taiwan § Department of Chemical Engineering, Rochester Institute of Technology, Rochester, New York 14623, United States ‡
ABSTRACT: The electrical double layer (EDL) profiles of electrical potential within charged slit-type pores is calculated in this work via a combination of grand canonical Monte Carlo (GCMC) simulations and electrodynamics concepts. The electrical potential distribution inside the pore is calculated with respect to field-free virtual bulk electrolyte solution implicit in the grand canonical ensemble. The GCMC simulations for slit-type pores performed in this work show that entropy effects lead to the exclusion of co-ions from the pore to allow densely packed counterions to efficiently neutralize surface charge. This phenomenon, in conjunction with the fact that electrical effects are long range with respect to the source charge, leads to a semioscillatory behavior of the potential profiles not predicted by classical theory. While bulk conditions in terms of ion concentrations are achieved in the center of the pore for most pore sizes studied in this work, electrical bulk conditions are not achieved except for the largest pores at the lowest surface charge density conditions. The electrical potential at the center of the pore is lower with respect to the absolute potential of the bulk solution in most cases.
I. INTRODUCTION
electrochemical reactions. Electrical double layers play a key role in applications related to energy storage and water purification. Moreover, such applications have shown dramatically improved performance when nanostructures are used.1−5 As the technology to manufacture nanodevices continues to improve,6 the study of EDL structure in such environments becomes increasingly relevant. Currently, there is no experimental method to directly probe EDL structure, so researchers must rely on mathematical theories or computer modeling to gain insights. The mathematical theory to describe EDL structure is known as the Gouy−Chapman− Stern (GCS) theory, based on the Poisson−Boltzmann
When a solid surface is immersed in a solution, the surface can become charged through a variety of natural thermodynamics processes or, in the case of an electrode, via the application of an external potential. The excess charge at the surface results in an accumulation of mobile ions in the solution near it. Relative to the charge on the surface, ions of opposite charge (counterions) build up near the surface to neutralize the surface charge, and ions of like charge (co-ions) are depleted in this region due to the effects of entropy and Coulomb repulsion. The resulting distribution of charges creates an excess electric potential at the interface. This potential drops to the Galvani potential of the bulk solution as the distance from the surface increases. The cumulative result is a phenomenon ubiquitous in electrochemistry known as the electrical double layer (EDL). The EDL governs natural phenomena like colloidal stability and electrokinetic phenomena; it also affects the outcome of © XXXX American Chemical Society
Special Issue: Emerging Investigators Received: November 30, 2017 Accepted: April 18, 2018
A
DOI: 10.1021/acs.jced.7b01048 J. Chem. Eng. Data XXXX, XXX, XXX−XXX
Journal of Chemical & Engineering Data
Article
equation.7 This classical theory assumes that the distribution of charge in the EDL is described by a Boltzmann distribution (with a slight modification proposed by Stern). This charge distribution is then used to solve the electrostatic Poisson equation to obtain both concentration and potential profiles. On the basis of how the energy in the Boltzmann distribution is defined, solutions to the GCS theory are implicit with respect to a zero potential7 in the bulk at a plane infinitely far away from the surface. More complex mathematical theories8−12 as well as numerous molecular modeling computer simulations13−18 have shown significant shortcomings in this theory, however, because it neglects the finite size of ions and accounts only for Coulombic interactions. Henderson and Boda give a comprehensive overview of mathematical theories and computer simulations for EDL and their relative successes and shortcomings.19 The present work focuses on the use of grand canonical Monte Carlo (GCMC) simulations, which can account for ion-size and ion− ion interactions, as a more robust means to study EDL structure. The grand canonical ensemble is considered ideal for simulating fluid inside a confined space that is in equilibrium with a much larger bulk, as is the case with nanopores, which are usually present in electrodes used in most electrochemical applications. There is a breadth of literature concerning the study of EDL concentration profiles using Monte Carlo simulations, both at a single charged interface and in confined systems.13,14,18,20−23 Contrary to GCS predictions, the observation of non-monotonic profiles and layering effects13,22,24 and attractive EDL interaction forces25,26 are noted. Other works explore the effects of mixtures of electrolytes. Such studies investigate the phenomenon known as ion selectivity, in which a specific ionic species may become immobilized in the EDL (i.e., electrosorbed) over other species. In these works, novel insights have been gained into phenomena such as charge inversion,27,28 competitive effects between ion size and ion valence,18,26,29−32 the effect of pore size on electrosorption capacitance,33−37 and the relationship between surface charge and ion selectivity.38,39 Many are continuing to explore the effects of nanopores and other systems on the scale of a few molecular radii on EDL structure using system models of varying complexity.40,41 There is comparatively less literature on obtaining full EDL potential profiles using molecular modeling. Most of the work in the open literature compute potential, ψ, via solving the onedimensional electrostatic Poisson equation d2ψ (x) 1 =− ρ (x ) ε0ε dx 2
ψ (x ) = −
∑ zieni(x) i
∫x
∞
(x′ − x)ρ(x′) dx′
(3)
This is typically integrated using numerical sums, which are computed either within the Monte Carlo code itself13 or using the output concentration profile to define the charge density as in eq 2.42,43 One work showed that this method is sensitive to statistical noise inherent in computer simulation results.44 Furthermore, the form of the solution in eq 3 is not suitable for confined systems, since it relies on boundary conditions that extend to infinity. Other studies have used the explicit Poisson equation integral solution ψ (x ) = −
1 ε0ε
x
⎡
x′
⎤
∫−∞ ⎢⎣∫−∞ ρ(x″) dx″⎥⎦ dx′ + C1x + C2 (4)
where C1 and C2 are constants of integration. However, one must determine the constants of integration to obtain a unique solution. The first constant of integration, C1, is defined via a known value of the electric field, typically at or beyond the walls of the simulation box. Defining the second integration constant, C2, is not so straightforward without using an arbitrary reference for potential. For example, one of the studies mentioned above sets the average potential in the center of the simulation box to zero in order to obtain absolute values of potential.45 Using this reference value can be ambiguous, particularly in small pores where the potential in the center of the pore can differ from values for the bulk solution. A third study solved the issue of arbitrary references by physically including the bulk solution in the simulation,45 via incorporating the charged walls and the bulk solution surrounding them inside the main simulation box. While this methodology is reflective of the system, it can be computationally very expensive. The goal of this work was to compute EDL potential profiles in confined systems using GCMC simulations. Potential profiles obtained from GCMC calculations were compared against classical GCS theory solutions and investigated the effects of confinement and surface charge on EDL potential profiles. In real systems, the potential applied to an electrode can be measured, but surface charge remains unknown. In molecular modeling, it is natural to prescribe a surface charge and see what the resulting EDL structure is. Hence, the relationship between surface charge and EDL structure is well documented. A method to compute EDL potential profiles from molecular modeling will therefore provide a link between surface charge, which defines EDL structure and interaction, and the experimental parameter applied potential.
(1)
This is accomplished using the concentration profiles (i.e., ion number density, n(x)), obtained from molecular modeling simulations. Here x describes the distance normal to the charged surface, ρ is the volumetric charge density, and ε and ε0 are the dielectric constant and permittivity of free space, respectively. The relationship between volumetric charge density and ion number density is given by
ρ (x ) =
1 ε0ε
II. METHODS II.A. Classical Theory Solution. The GCS theory can be used to obtain classical theory solutions for concentration profiles and potential profiles for the EDLs in slit-type pores. These solutions serve as a reference to which molecular modeling results can be compared and contrasted. The GCS theory assumes the structure of the charge distribution in the EDL. Gouy and Chapman capture the effects of mobile molecules by assuming the ion number density drops off exponentially as the distance from a charged wall increases according to the Boltzmann distribution. Stern modified this slightly, proposing that, if finite sized ions are modeled as hard spheres with point charges at their centers, they will have a distance of closest approach. Thus, no charge can exist less than
(2)
where z is the ion valence, e is the elementary charge, and subscript i represents individual ionic species. Many of the studies mentioned above use the one-dimensional convolution integral solution to eq 1 to obtain potential profiles. B
DOI: 10.1021/acs.jced.7b01048 J. Chem. Eng. Data XXXX, XXX, XXX−XXX
Journal of Chemical & Engineering Data
Article
deletion of an electroneutral pair of ions. The decision to pursue one of the moves over the other two was taken via the generation of a random number. Each type of move had a probability to be attempted equal to 1/3. Trial simulations were performed in order to determine the total number of cycles required for equilibration and averaging period for each case of surface charge density and pore diameter. During the equilibration period, the length of the random displacement of ions was evaluated every 10 cycles and was adjusted in order to keep the ratio of accepted moves between 25 and 30%. Parameter B, used in the partition function to evaluate insertion/deletion of ions, is a function of the excess chemical potential of the bulk electrolyte phase; i.e., it is a function of the activity coefficient and concentration of ions in bulk solution at equilibrium with the pore.28,29 Following the same practice of previous work, the activity coefficient and parameter B were determined from independent GCMC simulations of the bulk phase via an iterative procedure.28,29 The GCMC simulations of the bulk electrolyte phase were performed prior to the simulations of the electrolyte inside pores. The parameter B for the electrolyte concentration of 1.003 M was determined to be equal to 0.8803, and used as an input parameter in subsequent simulations. The temperature was equal to 298.15 K in all simulations. All simulations were performed on a Dell Precision T5610 workstation with an Intel Xeon Processor E5-2650 v2. More details on features and different algorithms of the custom-built Monte Carlo code are described in detail in previous works.18,26,28,29,31 The simulation box contains two impenetrable square walls of side length L at x = 0 and x = w, and periodic boundaries on the other two directions. A schematic of this system is shown in Figure 1. Table 1 outlines the system parameters explored in this
one molecular radius from a charged wall. Combining these two descriptions of the charge distribution and using eq 2, the volumetric charge density in a slit-type pore according to the GCS theory is written as ⎧ ⎪0 ⎪ ⎪ ⎛ zieψ ⎞ ⎪ 0 ⎟ ρ(x) = ⎨∑ zieni exp⎝⎜ kT ⎠ ⎪ i ⎪ ⎪0 ⎪ ⎩
0≤x
