Solution Structure Effects on the Properties of Electric Double Layers

Oct 24, 2018 - The structure of electrolyte solutions in electric double layers is primarily determined by the solvent, despite the fact that it is us...
0 downloads 0 Views 8MB Size
Subscriber access provided by Kaohsiung Medical University

Invited Feature Article

Solution Structure Effects on the Properties of Electric Double Layers with Surface Charge Regulation Assessed by Density Functional Theory Frank B. Van Swol, and Dimiter N. Petsev Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.8b02453 • Publication Date (Web): 24 Oct 2018 Downloaded from http://pubs.acs.org on October 28, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Solution Structure Effects on the Properties of Electric Double Layers with Surface Charge Regulation Assessed by Density Functional Theory Frank B. van Swol and Dimiter N. Petsev∗ Center for Micro-Engineered Materials and Department of Chemical and Biological Engineering, University of New Mexico, Albuquerque, NM 87131 E-mail: [email protected]

Abstract

1

The structure of electrolyte solutions in electric double layers is primarily determined by the solvent, despite the fact that it is usually neutral and not a subject to Coulombic interactions. The amount of solvent molecules in a typical electrolyte solution may be significantly greater that that of the ions. Hence, the charged species compete for space with a much larger number of neutral molecules, which has a strong effect on the density distributions near charged surfaces. Even for very dilute electrolyte solutions, the density profiles resemble liquid-like structure, which is entirely due to the presence of the dense solvent. Our work demonstrates that the solvent structural effect strongly couples to the surface chemistry, which governs the charge and potential. We argue that a comprehensive statistical-mechanical approach, such as the classical density functional theory that explicitly includes all solution species, in combination with a surface charge regulation condition at the interface, provides an excellent approach for describing charged interfaces. It allows for revealing important physical features and includes non-Coulombic contributions such as ionic and surface solvation.

Interfaces between electrolyte solutions and other substances occur in many systems of fundamental and practical importance. Examples include surface and colloid science, electrochemistry and corrosion, and soft materials and biomaterials. These interfaces are subject to complex physical and chemical interactions with the fluid, such as chemical reactions between surface terminal groups and certain solution components, physical adsorption, or desorption of charges at the interface. 1,2 All these interactions may lead to the accumulation of charge at the interface. As a result, the mobile charges in the vicinity redistribute in accordance with the electrostatic field generated by the surface. Hence, the effect propagates in the neighboring phases and is characterized by a spatial charge and electrostatic potential distributions. The whole system, consisting of charged interface and locally distributed mobile charges and potential is often referred to as the electric double layer (EDL). The name has historical roots, reflecting the initial notion that the charged interface with the electrolyte solution can be represented by a simple capacitor, where the surface charges and the solution counter-charges are placed on two welldefined planes in space. 3 This is certainly not an accurate representation of the physical sit-

Introduction

ACS Paragon Plus Environment

1

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

cently published by Trefalt et al. 2 The PB method suffers from a major deficiency. It assumes that the electrolyte solution is a structureless continuum and all effect, stemming from the actual size and various interaction between the ions and the solvent molecules are lost. This is particularly inaccurate at small (nanometer) length scales. Hence, very important fundamental and applied problems related to the properties of electrolyte solutions in confined spaces very near the charged interface cannot be adequately analyzed. The PB approach has been extended to account for the finite ion size and steric interactions near the charged interface, 17–19 but this approach remains still quite limited. A more accurate alternative is offered by modern statistical mechanics. However, its methods can be applied at different levels of detail. A popular model (sometimes referred to as “primitive”) takes into account the ion sizes and interactions, but continues to consider the solvent as structureless. For low to moderate electrolyte concentrations, the primitive model was shown to be equivalent to the modified PB approach discussed above. 20 A better model (known as “civilized”) should include the explicit structural contribution due to the solvent molecules and all interactions associated with them. 21–31 Note that the ”civilized” model would never approach the modified PB limit, with reducing the ion concentration, because of the structural effect of the solvent. The solvent density is much greater than the densities of all ionic species in the electrolyte solution. Hence, it dominates the overall liquid structure. Recently we examined the interplay between the local solution structure near the EDL charged interface, and the charge regulation chemical equilibrium that governs the surface charge using classical density functional theory (cDFT). 28–31 This approach is well-suited to account for a wide range of molecular and ionic interactions in the fluid, as well as with the adjacent substrate. In addition, cDFT has been extensively tested against computer simulations and it was demonstrated that the resultant density distributions from both approaches are in excellent agreement. 32 Our

uation, as all mobile charges (e.g, ions in the solution) are subject to Brownian motion at finite temperatures. An alternative model has been suggested, 4–6 based on Maxwell theory of electrodynamics, relating the potential Ψ(r) to the charge density distribution ρe (r) ∇2 Ψ = −

ρe , εε0

Page 2 of 24

(1)

where ε0 is the dielectric constant of vacuum, and ε is the relative dielectric permittivity of the solvent. Assuming an equilibrium, Boltzmann distribution of all ionic species   X −qi Ψ 0 . (2) ρe = ρi qi exp kB T i and inserting it in Eq. (1) leads to the very popular Poisson-Boltzmann (PB) model of the EDL. qi is the ion charge, and kB T is the thermal energy. The linearized version of the PB equation was used by Debye and Huckel to develop their celebrated theory of electrolyte solutions. 7 The PB model is based on the second order differential equation (1), and as such requires boundary conditions. The simplest examples are the Dirichlet (fixed potential), Neumann (fixed charge) conditions, or a combination of the two. This approach has the great advantage to be mathematically easy, and despite the nonlinearity of the PB model [see Eqs. (1) and (2)], allows to obtain analytical solutions for simple geometries. 4–6 However, a close inspection of the charging mechanisms of the interface with the electrolyte solution indicates that neither the Dirichlet nor the Neumann boundary conditions are justified. The surface charge is due to a thermodynamic equilibrium between the surface reactive groups and the solution components. A boundary condition that correctly accounts for the physical situation was first proposed by Ninham and Parsegian 8 and is often referred to as charge regulation. The surface charge regulation reflects the mechanism of charge formation and therefore can distinguish between different realistic cases. 9–16 An excellent overview of the PB analysis of EDLs using charge regulation boundary condition was re-

ACS Paragon Plus Environment

2

Page 3 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

cDFT analysis indicates that the solution structure is dominated by the presence of the neutral solvent and it strongly couples to the surface charge regulation chemistry. Hence, it has a significant impact on the EDL electrostatic properties. This is the main focus of the present study. A rigorous statistical mechanical approach allows to better interpret established quantities in electrochemistry, and surface and colloid science like the Helmholtz planes 3 and the Stern layer, 33 and offers a pathway for quantitative modeling of charged interfaces that involve electrolyte solutions.

2

energy, id

F [{ρi (z)}] = 2π kB T

   dz ρi (R, z) ln λ3i ρi (R, z) − 1 , (4)

p where λ = h2 /(2πmi kB T ) is the thermal de Broglie wavelength, h is Planck’s constant, mi is the mass of species “i”, and ρi (z) is the local density of component “i” along the z coordinate normal to the wall. The radial coordinate R runs parallel to the EDL interface with the substrate. The excess free energy consists of hard-sphere and long-range parts. The hard-sphere contribution is based on the derivation of Rosenfeld 36,37 and reads ex FHS [{ρi (z)}] = Z Z 2π kB T RdR dz ΦHS {nα (R, z)} . (5)

The grand thermodynamic potential

ΦHS {nα (R, z)} is the hard-sphere rediced free energy and nα (R, z 0 ) = nα (r) (r being the position vector) is the weighted local density.

EDLs involve multicomponent solutions, containing charged (ionic) and neutral (solvent) species. These species are subjected to the effects of the field exerted by the interface with the substrate. The field has electrostatic and van der Waals components, which lead to a variety of interactions with the fluid. In addition, the solution components (M ) in the fluid interact with each other. All of these interactions are easily incorporated in a cDFT model. The electrolyte solution is described in terms of a grand thermodynamic potential functional, which for single flat EDL reads, 27–31,34,35

ΦHS {nα (r)} = −n0 log(1 − n3 ) n1 n2 − n1 · n2 n32 − 3n2 n2 · n2 . (6) + + 1 − n3 24(1 − n3 )2 Remarkably, the functional form of the reduced free energy ΦHS {nα (R, z)} is independent on the number components M . It is the weighted densities nα (r) exhibit such a dependence and are defined by nα (r) =

Ω [{ρi (z)}] = id

RdR×

i=1

Z

Density functional theory of electric double layers with surface charge regulation

2.1

M Z X

M Z X

d3 r0 ρi (r0 )ωαi (r − r0 )),

(7)

i=1 ex FHS [{ρi (z)}]

ex [{ρi (z)}]+ Flong

F [{ρi (z)}] + + Z Z M X   2π RdR dz ρi (R, z) Viext (R, z) − µi .

and the weighting functions ωαi (r) are

i=1

(3) The first term on the right hand side of Eq. (3) corresponds to the ideal contribution to the free

ω3i (r) = Θ(Ri − r), ω2i (r) = δ(Ri − r), ω2i (r) i ω2i (r) i ω1 (r) = , ω (r) = , 4πRi 0 4πRi2 r ω i (r) ω i2 (r) = Θ(Ri − r), ω i1 (r) = 2 . (8) r 4πri

ACS Paragon Plus Environment

3

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The long-range contribution to the free energy functional is ex Flong [{ρi (z)}]

Z

ameter of component “i ”, and rij (R, z) is the distance between species “i ” and “j ”. All nonCoulombic interactions of a molecule (or ion) of type “i ” with the interface are assumed to be of the “hard wall” type, i.e. only the excluded volume effects are taken into account. The only exception is our analysis of the solvent-wall interactions presented in Section 3.3 below, where the LJ (9-3) potential 40 "    3 # 9 ds 2 ds , z > ds /2 − ΦLJ (z) = s−w 15 z z (14)

Z M M Z π XX RdR dz× = 2 i=1 j=1

dz 0 ρi (R, z)ρj (R, z 0 ) ΦLR (R, |z − z 0 |), (9)

where ΦLR (R, |z − z 0 |) is the long-range interactions contribution of the reduced free energy. Eq. (9) indicates that the long range interactions are accounted for in a mean-field limit. 38 The last term in Eq. (3) is the Lagrangian constraint, which accounts for the external fields Viext and fixed chemical potentials µi for all species. The densities of all components ρi (z) are found by minimizing the functional, 39 δΩ [{ρi (z)}] = 0. δρi (z)

is used in the analysis. The electrostatic interaction of the ions with the interface is Φel (z) =

ρe (z) =

qi ρi (z).

(11)

i=1

This approach includes the contributions of the interactions between the solutions species, which are present in the terms ex FHS [{ρi (z)}], 36,37 which accounts for the exex cluded volume effects, and Flong [{ρi (z)}], which captures long-range interactions. In our model, these interactions are of the Lennard-Jones (LJ) 40 "   6 # 12 dij dij ΦLJ (rij ) = 4ij − , rij > dij rij rij (12) and Coulombic Φel (rij ) =

qi qj e 2 , rij > dij , 4πεε0 rij

qi σz , z > di /2. 2εε0

(15)

Using a solution model based on a combination of Coulombic and LJ interactions is clearly an approximation when polar or polarizable solvents (i.e., water) are involved. Nevertheless, it captures important features of electrolyte solutions, such as attractive and repulsive interactions between all solution components. 28 This approximation works reasonably well for low to moderate ionic concentrations. 22,29,41–43 In such cases, the average distance between the charged species is sufficiently large so that any ioninduced local perturbation of the solvent molecular structure relaxes to the bulk state before reaching another ion. 44 For example in our case, the overall ionic concentration corresponds to 0.01 M. Assuming that the solvent molarity is 55.5 M, the average distance between two ions is ∼ 30 solvent molecular diameters. Even close to the charged surface, where the ion densities are greater (see Figure 1 below), the distance separating two ions is ∼ 10 molecular diameters. As Carnie and Chan showed, the local polarization perturbation decay to the bulk value over distances of the order of 3-4 molecular diameters. The overall free energy may also be reduced due to local polarization response of the solvent near the charged interface. We can mimic that effect by varying the solvent-wall interaction. While empirical and therefore not ideal, this approach has the advantage of being

(10)

The charge density distribution in the EDL is derived by summing over the individual charged species, taking into account the charge numbers with their sign M X

Page 4 of 24

(13)

type 27–31 where dij = (di + dj )/2, di is the di-

ACS Paragon Plus Environment

4

Page 5 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

The surface charge is 10

considerably simpler and computationally efficient. The Coulombic interactions are included in ex [{ρi (z)}] [see Eq. (9)], and are the term Flong therefore treated within a mean-field approximation. cDFT allows to include higher order correlations that are beyond the electrostatic mean-field approximation, but they lead to noticeable contributions at much higher ionic concentrations. 45

2.2

σ = eΓ

2

ρsAH + ρsAH+ + ρsA− 2

δ sinh[e(ΨN − Ψs )/kB T ] . (19) 1 + δ cosh[e(ΨN − Ψs )/kB T ] p where δ = 2 K− /K+ , Γ is the surface concentration of ionizable groups, and Ψs is the surface potential. The latter is directly related to the overall potential distribution in the EDL. Hence, even ions that are chemically inert and are not involved in the chemical reactions (16) affect the surface charge. This is due to simple charge screening, or more subtly, by affecting the chemical potential and therefore the local concentration of potential determining ions. The quantity ΨN = [ln(10)kB T /e] (pI − pHb ) is the Nernst potential where pI = (pK− +pK+ )/2 is the surface isoelectric point and pHb is the pH value in the bulk solution, far from the charged interface. The surface charge depends on the surface potential and the Nernst potential, and on the equilibrium constants K+ and K− . The Nernst potential depends on the pH in the solution bulk, and on the isoelectric point, which in turn depends only on K+ and K− . Relationship (19) corresponds to the charge at the physical interface between the substrate and the electrolyte solution, also known as the inner Helmholtz plane (IHP). 3 The charge of the Stern layer is located at the outer Helmholtz plane (OHP). 3,33 The actual physical situation could be more complicated as we discuss in Section 3.3. Still, for clarity we refer to the adjacent surface layer as the Stern layer when discussing some of the results below. Finding the density distributions of all ionic species in the solution allows for determining the overall charge density distribution more accurately than using the approximate Eq. (2). The charge density can be then introduced in Eq. (1) to find the electrostatic potential distribution in the solution and at the surface. The surface charge is also equal (and opposite in sign) = eΓ

Charge regulation and charge formation at the electric double layer interface

The electrostatic properties (charge and potential) at interface of the EDL is determined by a charge regulation boundary condition. 8,9 This condition derives the surface charge from the specific chemical equilibrium between the surface reactive groups and certain ionic species in the solutions – the potential determining ions (PDIs). This representation is much more realistic than fixing the surface charge or potential as external parameters. We use the model of Chan et al., 10 which considers the surface reaction equilibria + AH+ 2 + BH AH + BH2 , pK+ = − log10 K+ (16)

AH + BH A− + BH+ 2 , pK− = − log10 K− . in conjunction with the bulk reaction 2BH B− + BH+ 2

ρsAH+ − ρsA−

(17)

where K+ and K− are the equilibrium constants, and AH is a surface chemical group that can attach or release a proton depending on local densities of PDIs species BH+ 2 . This model is relatively simple, but it still allows to examine subtle effects such as pH variation around the surface isoelectric point. The densities of all ions in the immediate vicinity of the charged interfaces are given by 28,29 R di /2 dz ρi (z) s . (18) ρi = 0 R di /2 dz 0

ACS Paragon Plus Environment

5

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to the net bulk charge, or σ=−

M Z X i=1

10.56 and we use that number to determine the size of the system, which we set to 100d. It is also instructive to compare this longranged electrostatic length to the characteristic length-scale associated with the structural decay, which was found to be ∼ 1.3d. 28,29 The boundary condition for the surface charge is given by Eq. (19). The densities, far away from the interface, are equal to their fixed bulk values. More details of the computational procedure are provided elsewhere. 28,29



qi ρi (z)dz

(20)

0

Since the solvent is neutral, it is formally treated as a component with zero charge (qs = 0). We consider a flat EDL, which means that all charge distributions are averaged over planes parallel to the interface. This is evident from the integration over the radial coordinate R in Eqs. (4), (5) and (9). The charge regulation model, described above, also results in a surface charge which is uniformly distributed over the interface. In that respect, the surface charge treatment is in agreement with the cDFT representation of the adjacent fluid. Eqs. (1), (19), (11), and (20) are then iteratively solved to obtain the surface charges and potentials, as well as the corresponding distributions in the fluid phase. Defining a chemical equilibrium at the EDL interface (e.g., Eq. (19)) is equivalent to adding an appropriate surface functional term to Eq. (3). 46 The grand thermodynamic potential [see Eq. (10)] is minimized using the free-ware Tramonto code. 47 The parameters governing the surface chemistry are Γd2 /e = 0.66, pK− = 6, and pK+ = −2, which are common to oxide surfaces (e.g., silica). 28,29,48 This choice of parameters assigns a negative charge to the surface for all values of the solution pH greater than the pI = 2. The diameters of all solution species d are the same (except in cases when it is indicated otherwise). In order to ensure a correct conversion between scaled densities and molarity, we select d = 0.288 nm. 28,29 Thus, 1 L of solvent contains 55.5 M. The Debye screening parameter, is defined as 7 κ=

X 1 ρ0 q 2 εε0 kB T i i i

Page 6 of 24

3

Solution Structure and Surface Charge Regulation

The solvent in an electrolyte solution is involved in various interactions with the solutes (ions) and with the surface where the charge regulation reactions occur. These interactions are not of the simple Coulombic type, but they are still very important. The solvent occupies volume, which is unavailable to the ions. It may also be involved in attractive interactions with itself, the ions and the surface. Therefore, it is not surprising that the solution structure is significantly affected by these non-Coulombic solvent effect. The structure then determines the local densities of all species. Among those species are the PDIs, which take part in the surface chemical equilibria.(16) Hence, despite the neutral nature of the solvent molecules, they participate in a range of complex interactions, and because of their great number, may have a dominant effect on the properties of the EDLs. Below we discuss a few typical examples that are related to ionic and surface solvation.

!1/2 .

(21)

3.1

The Debye parameter has the dimension of inverse characteristic length that i relevant to the asymptotic charge and potential decay far from the interface. In the present study (κd)−1 =

Ion density distributions in the electric double layer

The ion density distributions in an EDL near the charged surface are shown in Figure 1. The plots are normalized by the corresponding bulk densities ρbi . The surface charge is governed by

ACS Paragon Plus Environment

6

Page 7 of 24

5

8

(a)

(z)/ bi

6

3 2

4 2

1 0

(b)

i

i

(z)/ bi

4

0

10

20

30

40

0

50

0

10

20

z/d

30

40

50

z/d

15

25

(c)

(d)

(z)/ bi

10

i

(z)/ bi

20

i

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

5

15 10 5

0

0

10

20

30

40

0

50

0

10

20

z/d

30

40

50

z/d

Figure 1: Ion density profiles ρi (z) near a charge regulating interface. ρbi is the bulk density of component i, far from the surface. The solid curves are for solvent diameter equal to that of the ions, d. The dashed curves are for solvent diameter equal to 0.01d. All blue curves correspond to the positive, non-PDIs, while the red curves – the negative ions. (a) pH = 4, I = 0.01 M; (b) pH = 4, I = 0.001 M; (c) pH = 5, I = 0.01 M; (d) pH = 5, I = 0.001 M. The remaining parameters are given in the text. Table 1: Surface potential and surface charge variations with the solution pH and electrolyte concentration. pH

I (M)

4 4 5 5

10−2 10−3 10−2 10−3

Structureless solvent eΨs /kB T σd2 /e −0.883 − 0.00264 −1.440 − 0.00147 −2.222 − 0.00766 −2.872 − 0.00369

Explicit eΨs /kB T −0.721 −1.239 −1.980 −2.634

solvent σd2 /e − 0.00215 − 0.00122 − 0.00666 − 0.00322

to a larger resultant surface charge. The variation of the bulk density of ions, or the ionic strength, primarily has an effect on the range of the electrostatic potential into the solution phase. However, since the ionic strength mod-

the regulation condition (19). We explore the combination of two pH values and two different ionic strengths. The pH variation affects the surface charge, i.e., a greater difference between the solution pH and the surface pI leads

ACS Paragon Plus Environment

7

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ifies the local potential distribution, it leads to slight variations in the surface charge [see Eq. (19)]. The selected surface parameters and pH values result in a negative surface charge. Hence, the positive ions (blue curves) are attracted to the interface and their density locally increases. The negative ions (red curves) are repelled and hence, their density near the charged interface is lower. The solid curves in Figure 1 show the density distributions obtained from the cDFT, explicitly taking into account the presence of the solvent molecules. The diameter of all species is the same. All interactions between the solution components include excluded volume, and LJ attraction with energy parameter ij = 1 kB T [see Eq. (12)]. In addition, all ions experience Coulombic interactions in accordance with Eq. (13). It is important to emphasize that the ionic density is rather low (gas-like), but the density profiles exhibit well defined (liquid-like) structural peaks. This structuring is due to the explicit presence of the solvent, which dominates the underlying liquid structure. The ions can only be placed in a void space that has not been occupied by a solvent molecule. This means that they effectively replicate the dense liquid structure. The electrostatic field exerted by the charged interface contributes to the slopes of the curves at larger distances. These slopes are overall negative for the positive ions and positive for the negative ions. The dashed curves illustrate the electrolyte behavior in the case where the neutral solvent is practically structureless (i.e.,the primitive model limit). The detailed interactions involving the solvent molecules are not ignored. At the same time, the ion densities are too low to produce any noticeable structural peaks near the charged surface. Another observation is that the curves with and without the solvent contributions do not perfectly match, even at distances where the structure is no longer present. This is due to the solventinduced structural effect on the surface charge regulation [see eq. 19]. The surface charges and potentials are slightly more negative when the solution does not exhibit any structuring

Page 8 of 24

as shown in Table 1. This finding demonstrates that ignoring the solvent-induced solution structure leads to overestimated results for magnitude of the surface charges and potentials. Using traditional boundary conditions such as constant charge or potential at the surface would also miss this effect. Hence, it is the coupling between solution structure and surface chemistry that reveals this detail. Figure 1a shows the ionic density profiles for pH = 4 and I = 0.01 M. The surface charge and potentials are low (see Table 1), and the Debye screening length is equal to κ−1 = 10.56d. For an ionic diameter of 0.288 nm, the dimensional screening corresponds to 3.04 nm. Reducing the ionic strength down to I = 0.001 M increases the screening length to κ−1 = 33.38d, or 9.61 nm. This is the reason for the more gradual decrease in the local density with distance (see Figure 1b). In addition, the background electrolyte concentration reduction also leads to an increase of the magnitude of the surface potential and a decrease of the magnitude of the surface charge (see Table 1). Changing the solution pH to 5 (while keeping I = 0.01 M) increases the magnitudes of both the surface potential and charge (see Figure 1c and Table 1). As a result, the local density of positive ions near the surface increases while that of the negative ions decreases. Finally, Figure 1d depicts the case where pH = 4 and I = 0.001 M. The surface potential magnitude increases further, while the surface charge magnitude drops in comparison with the higher electrolyte concentration case shown in Figure 1c. The local density of positive ions near the surface is even higher, and the one for the negative ions is lower. The propagation of the charge and potential into the fluid is greater because of the larger screening length. An important conclusion stemming from the above analysis is that the explicit account for the neutral solvent effects leads to a proper description of the solution structure in the EDL. This structure is important for the correct computation of the surface potential and charge in the presence of surface charge regulation [see Eq. (19)]. The surface potential and charge then determine the properties of the whole EDL

ACS Paragon Plus Environment

8

Page 9 of 24

vation energy (i.e., the LJ parameter s−PDI ). As the PDIs become more solvated, their density near the interface with the substrate decreases and the peaks move farther into the solution. The reduction of PDI concentration shifts the chemical equilibrium (17) towards a more negative surface [see also Eq. (19), and Figs 3a and 3b]. This leads to a strong attraction of positive charge in the fluid in the immediate vicinity of the negatively charged interface, as shown by the peaks in the fluid charge density distribution in Figure 3c. The charge distribution in the EDL is naturally accompanied by an increase in the potential distribution (see Figure 3d). The surface charge and potential exhibit a particularly abrupt change around s−PDI /kB T ∼ 3.4. This change correlates with a similarly abrupt shift in the PDI density distributions shown in Figure 2. For the selected parameters listed above, σe/d2 = 0.664 corresponds to a fully ionized surface (see Figure 3a). The surface potential magnitudes that correspond to a fully ionized surface are nonphysically high, which implies that such a scenario is unlikely. Still, it clearly shows how the surface responds to the extreme solvation of the PDIs through the surface charge regulation mechanism. The electrostatic potential distributions in the EDL (see Figure 3d) exhibit a linear change between z/d = 0 and 0.5. The plane z/d = 0.5 corresponds to the location of the center of charges in the first layer of fluid next to the charged interface (see Figure 3c). This plane defines the location of the Stern layer. 33 More discussion of the the Stern layer in the framework of cDFT applied to charge regulating EDLs is provided below (see also Refs. 29–31 ).

50

/ bPDI

40 30

PDI

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

20 5

10

4 3

0

2 0

2

4

z/d

6

8

10

1

/k T

s-PDI B

Figure 2: PDI density profiles for varying solvation LJ energy s−PDI /kB T . including the far regions where structure has already disappeared. The solution structure varies with the distance from the surface as ∼ exp(−βz) with characteristic length β −1 = 1.3d. 28,29 Taking into account the structuring near the surface allows for determining the composition of the Stern layer 30,31 and the location of the electrokinetic shear plane. 29

3.2

Effects due to ionic solvation

An important advantage of solvent-explicit molecular models of EDLs is that they can naturally incorporate the variety of non-Coulombic interactions that involve the neutral solvent molecules and the ionic species in the solution. These effects can be extremely important for the properties of electrolyte solutions and EDLs. 49–51 The solvation interactions can be relatively simply taken into account in the framework of a LJ solvent. Clearly this is an approximation for the case of polar solvent such as water, but it is a reasonable one, provided that the ionic concentration is low to moderate. 28,29 We have recently studied the effect of weak and average solvent-solvent and solvent-ion interactions, and found that they have a significant effect properties of charged interfaces. 30 Below we offer new results that extend the analysis to the strong ionic solvation limit. 3.2.1

3.2.2

Solvation of the positive ions

The presence of background electrolyte (consisting of ions that are not directly involved in chemical, potential determining interaction with the interface) may still have a strong effect on the properties of EDLs. These ions may still be involved in physical (Coulombic, or LJ) interactions with the surface as well

Solvation of the potential determining ions

Figure 2 shows the variation of local PDI density with distance and the strength of their sol-

ACS Paragon Plus Environment

9

Langmuir 0

-0.3

(b)

(a) -5

e s/kBT

d2/e

-0.4

-0.5

-0.6

-0.7

-10 -15 -20 -25

1

2

3

4

/k T s-PDI B

5

1

2

3

4

/k T

5

s-PDI B

0 (d) -5 (c)

0.4

0.2

6

0.1

4

0.0 0

e

d3/e

e

d 3/e

e /kBT

0.3

8

2

4

z/d

5 4 3 2 s-PDI /kB T 1

-10 -15 -20

2 5 4

0

3 2 0

2

4

z/d

6

8

10

/k T

-25

s-PDI B

1

0

2

4

6

8

10

z/d

Figure 3: Effect of the PDI solvation on the (a) surface charge; (b) surface potential; (c) total charge distribution in fluid near the charged surface (the inset provides an enlarged view); (d) potential distributions in the fluid (the potentials become more negative with s−PDI – see Figure 3b). The linear potential regions between z/d = 0 and 0.5 correspond to the distance between the charged interface and the centers of fluid molecules in the first adjacent layer [see also the first peak in the charge distribution in (c)]. Adapted from J.Chem.Phys. 2017, 147, 214704; DOI: 10.1063/1.5005060. as with all species in the fluid including the PDIs. The positive ions carry the same charge as the PDIs, and therefore, they are attracted by the negative surface. Hence, they compete for space in the Stern layer, which affects the surface charge by shifting the equilibrium (17). The density profiles for the positive ions near the charged interface and for different solvation energies s−Pos , shown in Figure 4. The positive ion density in the Stern layer decreases with the increase of the solvation energy. For s−Pos & 2kB T , the positive ions are practically absent from the layer in the immediate vicinity of the charged surface. Instead, they form layers of ions farther into the fluid as indicated by the structural peaks. This observation is in qualitative agreement with X-Ray measurements of effect of hydration on the local ionic densities near muscovite-water interfaces. 52 The solvation of the positive ions has an effect on the charge and potential at the surface and

5

/ bPos

4 3

Pos

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 24

2 5

1

4 3

0

2 0

2

4

z/d

6

8

10

/k T

s-Pos B

1

Figure 4: Positive ions density profiles for varying LJ solvation energy s−Pos /kB T

ACS Paragon Plus Environment

10

Page 11 of 24

-1.9

-0.70

(b)

-2.0

-0.75

e s/k BT

d2/e 103

(a)

-2.1

-0.80

-2.2

-0.85 1

2

3

4

5

1

2

3

/k T

4

5

/k T

s-Pos B

s-Pos B

-0.2

(c)

e /kBT

d3/e 104

(d)

-0.4

4

e

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

2

-0.6

-0.8

0

5 4 3

-2

2 0

2

4

6

z/d

8

10

1

s-Pos/kBT

-1.0

0

2

4

z/d

6

8

10

Figure 5: Effect of the positive ions solvation on the (a) surface charge; (b) surface potential; (c) total charge distribution in fluid near the charged surface; (d) potential distributions in the fluid (the potentials become more negative with s−Pos – see Figure 5b). The linear potential regions between z/d = 0 and 0.5 correspond to the distance between the charged interface and the centers of fluid molecules in the first adjacent layer [see also the first peak in the charge distribution in (c)]. Adapted from J.Chem.Phys. 2017, 147, 214704; DOI: 10.1063/1.5005060.

ACS Paragon Plus Environment

11

Langmuir

in the fluid (see Figure 5). For s−Pos = 1kB T the interaction is the same as for the remaining fluid components, including the PDIs. This means that both the positive ions and the PDIs interact with the wall and the solvent exactly the same. Both tend to occupy the Stern layer, and because of their identical charge, they experience mutual repulsion. Hence, the positive ions create a barrier for the PDIs, which is both steric and electrostatic, and reduce their concentration in the layer next to the surface (see also Refs. 29,30 ). This leads to a highly negative charge at the surface. As the positive ions become more solvated, they are removed from the layer that is in immediate contact with the surface. Hence, the PDIs have better access to the reactive groups and reduce the negative charge by means of the charge regulation reaction (17), as shown in Figure 5a. However, as the ions are removed from the region near the

3.2.3

1.0

/ bNeg

Solvation of the negative ions

Figure 6 shows the density profile for the negative ions in the EDL for varying solvation energy s−Neg . Similar to the cases for the PDIs and the positive non-PDIs above, the local negative ion density shifts away from the surface into the fluid as the solvation energy increases. The overall effect of the negative ions solvation on the surface charge and potential, as well as on their distributions in the fluid is relatively weak (see Figure 7). The surface charge (Figure 7a) becomes slightly more negative as the negative ions become more solvated. At the same time, the electrostatic surface potential (see Figure 7b) in the fluid becomes (again slightly) less negative. Consequently, positive ions move in close to the surface which results in a net positive fluid charge in the Stern layer and beyond (see Figure 7c). The potential distribution in the EDL (see Figure 7d) follows the trend determined by the surface values, and the screening due to the dissolved ions. The linear portions of the curves between z/d = 0 and 0.5 are again due to the potential change inside the Stern layer.

1.2

Neg

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 24

0.8 0.6 0.4 5

0.2

4 3

0.0

2 0

2

4

z/d

6

8

10

1

/k T

3.3

s-Neg B

Effects due to solvent-surface interactions

The interface between a charged wall and the fluid in an EDL can be considered as a system component similarly to those in the fluid. It interacts with all solution species, including the solvent. Here we describe the effect of longranged solvent interaction with the wall in the framework of the LJ model [see Eq. (15)]. Below we examine the effect of the solvent-wall interactions of the LJ type, while the ions experience only Coulombic and excluded volume forces. Analysis of the effects of the LJ ionwall attractive interactions and their competition with the solvent on the properties of EDLs is presented elsewhere. 31 Figure 8 shows the density profiles for the solvent near the charged interface. The different curves are for varying attractive strength with the wall. The first peak in the solvent density exhibits a noticeable increase with LJ energy

Figure 6: Negative ions density profiles for varying LJ solvation energy s−Neg /kB T charged wall, the potential becomes more negative (see Figure 5b). The charge distribution in the fluid is shown in Figure 5c. As the solvation of the positive ions increases (with s−Pos ), the Stern layer becomes more populated with negative ions. Farther away from the surface, however, the fluid charges exhibits positive peaks. Hence, the solvation effects lead to strong structuring of both the liquid and the charge near the surface. The electrostatic potential distribution in the EDL again exhibits a linear portion between z/d = 0 and 0.5, and monotonically increases beyond the Stern layer (see Figure 5d).

ACS Paragon Plus Environment

12

Page 13 of 24

-2.18

-0.68

(b)

(a) -0.69

e s/k BT

d2/e 103

-2.19

-2.20

-2.21

-0.70 -0.71 -0.72

-2.22 1

2

3

4

-0.73

5

1

2

3

/k T

4

5

/k T

s-Neg B

s-Neg B

0.0

(d)

(c)

-0.5

e /kBT

d3/e 104

8 6 4

-1.0

e

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

2

5 4 3

0

2 0

2

4

z/d

6

8

10

/k T

s-Neg B

-1.5

1

0

2

4

z/d

6

8

10

Figure 7: Effect of the negative ions solvation on the (a) surface charge; (b) surface potential; (c) total charge distribution in fluid near the charged surface; (d) potential distributions in the fluid (the potentials become less negative with s−Neg – see Figure 7b). The linear potential regions between z/d = 0 and 0.5 correspond to the distance between the charged interface and the centers of fluid molecules in the first adjacent layer [see also the first peak in the charge distribution in (c)]. Adapted from J.Chem.Phys. 2017, 147, 214704; DOI: 10.1063/1.5005060.

ACS Paragon Plus Environment

13

Langmuir

charge peaks in Figure 9c are also very close to linear. These linear portions of the potential distributions correspond to the changes between the layer ions near the surface.

10

/ bs

8 6

s

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 24

4 2 0 0

2

4

z/d

6

8

10

5 4 3 2 /k T 1 s-w B 0

3.4

Figure 8: Solvent density profiles for varying LJ solvent-charged wall interactions s−w /kB T

Molecular interpretation of the Stern layer and Helmholtz planes

It is useful to relate the results obtained by the cDFT approach to established concepts from traditional electrochemistry, and colloid and interface science. According to Stern, 33 the first layer that is adjacent to the charged surface has a special significance and merits a separate analysis. This physical representation is very similar to the first Helmholtz model of the EDL. 3 Stern, however, considered the region outside of the adjacent layer to be described by the Gouy-Chapman model. 4–6 Hence, the Stern layer is modeled as an adsorption, while the remaining EDL is analyzed by means of the Poisson-Boltzmann equation [i.e., the combination of Eqs. (1) and (2)]. Grahame 53 refined the concept by defining a plane where the physically adsorbed or covalently bound ions are located, often referred to as the inner Helmholtz plane. Fully solvated ions form the next, outer Helmholtz plane. Clearly, these ideas are oversimplified. While correctly assuming that the ion dimensions are important, they are not taken into account in the diffuse part of the EDL. In addition, factors such as surface roughness or ion-size variations do not allow the placement of the ions on well-defined planes. Finally, the effect of the solvent is underestimated, since it plays a dominant role in the local structuring near the charged interface. A statistical mechanical model with molecular resolution, like the cDFT discussed and applied in this paper, does not require an ad hoc introduction of particular structure such as layers and planes. The molecular organization of the solution in the EDL comes out naturally as a result of minimization of the thermodynamic potential [see Eq. (10)]. Still, a layered structure of all the species (charged ions and neutral solvent) emerges, undoubtedly enhanced by the

parameter s−w . The second, third, and the remaining distant peaks also become more pronounced, which indicates that the fluid becomes more structured (i.e., layered) near the wall. This has an effect on the local density distribution of all ionic species, which in turn affects the charges and potentials both at the surface and in the fluid as shown in Figure 9. The surface charge (see Figure 9a) increases with the local density of the solvent molecules, which displace, and effectively dilute the concentration of PDIs. Thus, the surface equilibrium (16) is shifted towards greater dissociation of the reactive groups. The surface potential (see Figure 9b) also becomes more negative with the LJ parameter s−w . The charge distribution in the fluid is shown in Figure 9c. Interestingly, it indicates some positive charge density in the adjacent layer to the wall, but a lot more in the next few layers. Hence, while the ions are mostly displaced near the surface, some remain there attracted by the Coulombic force exerted by the oppositely charged interface. As the solvent density drops with the distance from the wall (see Figure 8), the (mostly positive) ions form structured layers, which are again due to the attractive and long-ranged Coulombic force (see Figure 9c). The potential distribution curves in the EDL and for different values of s−w are plotted in Figure 9d. The potentials are negative and exhibit a well-defined linear region between z/d = 0 and 0.5. A close inspection of the potential curves reveals that the changes between the regions corresponding to the position of the main

ACS Paragon Plus Environment

14

Page 15 of 24

Langmuir

0.0

0

(a)

(b) -1 -2

e s/kBT

d2/e 102

-0.5 -1.0 -1.5

-3 -4

-2.0 -2.5

-5 0

1

2

/k T s-w B

3

4

-6

5

0

1

2

/k T

3

4

5

s-w B

0 -1 (c)

e /kBT

/ b 103

(d)

-2

6 6

-3 -4

4

e

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2 0 0

2

4

z/d

6

8

10

-5

5 4 3 2 /k T 1 s-w B 0

-6

0

2

4

6

8

10

z/d

Figure 9: Effect of the solvent wall-surface interactions on the (a) surface charge; (b) surface potential; (c) total charge distribution in the fluid near the charged surface; (d) potential distribution in the fluid (the potentials become more negative with s−w – see Figure 9b). The linear potential regions between z/d = 0 and 0.5 correspond to the distance between the charged interface and the centers of fluid molecules in the first adjacent layer [see also the first peak in the charge distribution in (c)]. Adapted from J.Chem.Phys. 2018, 148, 044702; DOI: 10.1063/1.5012090. fact that we also assumed a molecularly smooth interface between substrate and solution, 54 as well as because all solution species are assumed to have the same size. Our analysis suggests that indeed there are certain important planes where species tend to locate. The first one is the plane where surface charge regulation chemical reaction occurs. The model we used above distinguishes this plane from the one where physical adsorption of ions or solvent molecules may occur. This is different from Grahame’s representation, where both the chemical and physical bonding occurs on the same plane. There are common realworld examples which support such a distinction. An oxide surface in contact with aqueous solution will regulate its charge by binding or releasing protons. However, these are usually carried around as hydronium ions, which are much larger. Therefore, the plane the hydronium ions reach as they touch the surface is not the same as the reaction plane, where the

chemical interaction involving the much smaller hydrogen ion occurs. In addition, the layered structure can persist for more than a single layer of ions and molecules. Hence, the realistic physical picture can be more complex that the simple Stern, or even Graham’s models, and is entirely governed by the particular type and strength of the interactions between all solutions species and with the charged surface (see also Refs. 30,31 ) In light of this discussion, it is useful to consider the relationship between the charge and potential distributions in the fluid. While the fluid charge distributions exhibit multiple peaks near the surface, the respective potential distributions are mostly monotonic. This is not surprising since Eq. (1) states that the charge distribution is related to the curvature of the potential function. Hence, as long as the charge distribution does not change sign, the curvature of the electrostatic potential distribution functions will not change either. Obviously, a con-

ACS Paragon Plus Environment

15

Langmuir

stant curvature requirement prevents the potential distribution from exhibiting the oscillatory behavior that is observed for the charge or the density profiles. There are examples, like in Figure 5c above, where the fluid charge density changes sign, and yet the corresponding potential curves, shown in Figure 5d, remain monotonic (more examples can be found in Refs. 30 and 31 ). Close inspection of the potential curves, however, shows that there are variations in the curvatures, but their magnitude is too small to be visible at this scale. This is also why the potential curves do not exhibit minima and maxima like the fluid charge distributions. In fact, differentiating twice the potential curves in Figure 5d recovers the charge density curves in Figure 5c. Another way to demonstrate the curvature variation in the potential distribution is by subtracting the linear contribution, which dominates for small z/d, and examine the nonlinear residual, ∆ΨN L defined by     ∂Ψ z . ∆ΨN L (z) = Ψ(z) − Ψs + ∂z z=0 (22)

0.5 and 1.35 and negative at larger distances. Hence, the curvature variations in the potential distribution function are in perfect agreement with the results for the fluid charge [see also Eq. (1)]. Note that Eq. (1) is as accurate as charge density distribution ρe on its right hand side. Combining Eq. (1) with Eq. (2) leads to the rather approximate Poisson-Boltzmann model. Having a more accurate distribution that accounts for size and molecular and ionic interactions leads to a much more realistic relationship between fluid charge and potential. Such a framework allows one to analyze the entire EDL as a single entity, instead of rather artificially breaking it into layers that are governed by different models, as suggested by Stern 33 and Grahame. 53

4

-3.0

NL B

/k T

102

-2.9

-3.1

-3.2 0.0

0.5

1.0

1.5

2.0

2.5

Conclusions

Electric double layers are complex and diverse systems which present many interesting fundamental physical and chemical problems. At the same time, they appear as important elements in many established and emerging practical applications. Therefore, it is not surprising that electric double layers have been extensively studied for more that a century. Still, quite often the theoretical models that are used are rather oversimplified. In a recent Feature Article, Trefal at al. 2 presented a strong case in favor of using a surface charge regulation condition at the substrate-electrolyte interface instead of setting a constant potential or charge as an external parameter. In the present paper, we argue that an account for the precise solution structure is extremely important. Furthermore, we demonstrate [see also Refs. 28–31 ] that this structure is largely determined by the neutral solvent, despite the fact that it does not participate in any Coulombic interactions. Typically, the number of solvent molecules greatly exceeds that of the ions, which is the reason for its strong effect. The charged ionic species are forced to assume locations in a structural matrix created by the dense solvent. That is why their density profiles near the surface resemble a liquid-like structure

e

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 24

3.0

z/d

Figure 10: Nonlinear residual of the potential distribution [see eq. (22)] in the EDL for s−Pos = 2kB T The quantity ΨN L is plotted in Figure 10. The plotted curve is for s−Pos = 2kB T , which is shown in Figure 5c, corresponds to a case where the fluid charge density exhibits a negative minimum between z/d = 0.5 and 1.35. For z/d > 1.35, the charge density is always positive. The plot in Figure 10 shows that the potential curvature is positive between z/d =

ACS Paragon Plus Environment

16

Page 17 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

while their actual concentration is gas-like. The solution structure near the charged interface affects the surface chemistry, and hence there is a strong coupling between the solventinduced effects and the charge regulation [for more details see Refs. 28 and 29 ]. The surface charge then determines the electrostatic force experienced by the ions in the electric double layer. The intermolecular interactions inside an electrolyte solution are complex. Traditionally, in an effort to improve our understanding of solution phenomena we break these interactions up into a variety of contributions, e.g., we speak of ionic, polar, dipolar, induced dipolar and hydrogen bonding interactions. It is understood, but not always stressed, that the division of these contributions is arbitrary. That is, experimentally there are no ways to consider physical quantities, such as for example the heat of dissolution of a salt, say, and assign the ionic and polarization contributions to that quantity. The situation is different for modeling approaches, such as MD and cDFT, because there we can specify a Hamiltonian and allow ourselves to perform thought experiments. Traditionally, we attempt to keep the design of a Hamiltonian as simple as possible in order to facilitate the understanding and physical intuition that is to be generated by the study of model systems. The celebrated van der Waals model for fluids and fluid mixtures is a good example in this context. It consists of a just two kinds of contributions: excluded volume (hard sphere interactions) and a mean field attractions. Despite its simplicity, the van der Waals model has been extremely effective in the description of, say, phase diagrams and interfacial phenomena. However, the study of critical phenomena demonstrated that the value of the critical exponents were not quite right. Our approach toward the electrolyte solutions proceeds in a similar vein. Clearly, we recognize that the inclusion of the solvent introduces a plethora of interactions, as indicated above. However, at this stage we have sought to keep the computational challenges manageable by simplifying the interactions. It is in this spirit that we consider representing the solvent

by a simple Lennard-Jones interaction, which therefore to be considered an effective representation of water. This, we believe, captures the most important aspects: excluded volume and attractive interactions. The ion-solvent interactions are generally of the order of the thermal energy. Poorly solvated ions have parameter that is less than KB T , while the parameter for strongly solvated ions is larger than kB T . This is systems-specific, and that is why we covered a range of values to examine how different cases would affect properties like charge and potential distributions. A model that explicitly includes the solvent at the molecular level allows one to examine the effect of solvation non-Coulombic interactions. These are very common and are responsible for the variations in the electric double layer properties containing ions with the same ionic charge, but different otherwise. 49–51 In this paper we demonstrated the power of cDFT to model the effects of the solvation interactions in combination with the surface chemical equilibrium on the charge and potential distributions at the surface and in the fluid. More examples are available elsewhere. 30,31 A detailed molecular model of the electric double layer offers a better physical interpretation of established concepts such as the Stern layer and the Helmholtz planes. An elaborate molecular theory demonstrates that these concepts, while sometimes useful, are usually oversimplified and the real physical picture is more complex. The charge layering has been known for a few decades and attempts were made to incorporate it into a theoretical model. Many of these models, however, ignore the effect of the solvent contribution to the solution structure. Some of these models are based on modifying the Poisson-Boltzmann model, 17–19,55,56 or employed integral equations approaches, while assuming that the solvent is a structureless continuum. 57–62 Other statistical-mechanical analyses did include the solvent molecular contributions, but did not account for the surface charge regulation. 21–26,63 The message of the present study is that the solvent structure and the surface chemistry are coupled, and taking one into account while ignoring the other may result in

ACS Paragon Plus Environment

17

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

missing crucially important properties of electric double layers.

Page 18 of 24

(8) Ninham, B. W.; Parsegian, V. A. Electrostatic Potential between Surface Bearing Ionizable Groups in Ionnic Equilibrium with Physiologic Saline Solution. J. Theor. Biol. 1971, 31, 405–428.

Acknowledgement This research was supported by the Air Force Academy under cooperative agreement FA7000-14-2-0017 through funding from the Corrosion Policy and Oversight office, and the United States Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering. The authors are indebted to Professor David G. Whitten for his encouragement and comments on an earlier version of the manuscript. We are also thankful to the UNM Center for Advanced Research Computing for providing support and access to computational resources.

(9) Ettelaie, R.; Buscall, R. Electrical double layer interactions for spherical charge regulating colloidal particles. Adv. Coll. Interf. Sci. 1995, 61, 131–160. (10) Chan, D. Y. C.; Healy, T. W.; White, L. R. Electrical Double Layer Interactions under Regulation by Surface Ionization Equilibria–Dissimilar Amphoteric Surfaces. J. Chem. Soc., Faraday Trans. I 1976, 72, 2844–2865.

References

(11) Chan, D. Y. C.; Healy, T. W.; Supasiti, T.; Usui, S. Electrical double layer interactions between dissimilar oxide surfaces with charge regulation and Stern– Grahame layers. J. Colloid Interface Sci. 2006, 296, 150–158.

(1) Dukhin, S. S. In Surface and Colloid Science; Matijevic, E., Ed.; Wiley Interscience: New York, 1974; Chapter 1, pp 1–49. (2) Trefalt, G.; Behrens, S. H.; Borkovec, M. Charge Regulation in the Electrical Double Layer: Ion Adsorption and Surface Interactions. Langmuir 2016, 32, 380–400.

(12) Behrens, S. H.; Grier, D. G. The charge of glass and silica surfaces. J. Chem. Phys. 2001, 115, 6716–6721. (13) Borkovec, M. Origin of 1-pK and 2-pK Models for Ionizable Water–Solid Interfaces. Langmuir 1997, 13, 2608–2613.

(3) Helmholtz, H. Ueber einige Gesetze der Vertheilung elektrischer Strme in krperlichen Leitern mit Anwendung auf die thierisch-elektrischen Versuche. Ann. Phys. Chem. 1853, 165, 211–233.

(14) Behrens, S. H.; Borkovec, M. Electric double layer interaction of ionizable surfaces: Charge regulation for arbitrary potentials. J. Chem. Phys. 1999, 111, 382–385.

(4) Gouy, G. Sur la constitution de la charge electrique ‘a la surface dun electrolyte. J. Physique 1910, 9, 457–468.

(15) Behrens, S. H.; Borkovec, M. Electrostatic Interaction of Colloidal Surfaces with Variable Charge. J. Phys. Chem. 1999, 103, 2918–2928.

(5) Gouy, G. Sur la fonction lectrocapillaire. Ann. Phys 1917, 7, 129–184. (6) Chapman, D. L. A contribution to the theory of electrocapillarity. Phil. Mag. 1913, 25, 475–481.

(16) Popa, I.; Sinha, P.; Finessi, M.; Maroni, P.; Papastavrou, G.; Borkovec, M. Importance of Charge Regulation in Attractive Double-Layer Forces between Dissimilar Surfaces. Phys. Rev. Lett. 2010, 104, 228301.

(7) Debye, P.; Huckel, E. Zur theorie der elektrolyte. I. Gefrierpunktserniedrigung und verwandte erscheinungen. Phys. Ztschr. 1923, 24, 185–206.

(17) Borukhov, I.; Andelman, D.; ; Orland, H. Steric effects in electrolytes: A modified

ACS Paragon Plus Environment

18

Page 19 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Poisson-boltzmann equation. Phys. Rev. Lett. 1997, 79, 435–438.

(27) Frink, L. J. D.; van Swol, F. Solvation Forces and Colloidal Stability A Combined Monte-Carlo and DensityFunctional Theory Approach. J. Chem. Phys. 1994, 100, 9106–9116.

(18) D. Ben-Yaakov, D. H., D. Andelman; Podgornik, R. Beyond standard PoissonBoltzmann theory: Ion-specific interactions in aqueous solutions. J. Phys.: Condens. Matter 2009, 21, 424106.

(28) Fleharty, M. E.; van Swol, F.; Petsev, D. N. Solvent Role in the Formation of Electric Double Layers with Surface Charge Regulation: A Bystander or a Key Participant? Phys. Rev. Lett. 2016, 116, 048301.

(19) Ben-Yaakov, D.; Andelman, D.; Podgornik, R.; Harries, D. Ion-specific hydration effects: Extending the PoissonBoltzmann theory. Curr. Opin. Colloid Interface Sci. 2011, 11, 542–550.

(29) Vangara, R.; Brown, D. C. R.; van Swol, F.; Petsev, D. N. Electrolyte solution structure and its effect on the properties of electric double layers with surface charge regulation. J. Colloid Interface Sci. 2017, 488, 180–189.

(20) Attard, P. Electrolytes and the Electric Double Layer. Adv. Chem. Phys. 1996, 92, 1–159. (21) Adelman, S. A.; Deutch, J. M. Exact solution of the mean spherical model for strong electrolytes in polar solvents. J. Chem. Phys. 1974, 60, 3935–3949.

(30) Vangara, R.; van Swol, F.; Petsev, D. N. Solvation effects on the potential and charge distributions in electric double layers. J. Chem. Phys. 2017, 147, 214704.

(22) Chan, D. Y. C.; Mitchell, D. J.; Ninham, B. W.; Pailthorpe, B. A. On the theory of dipolar fluids and iondipole mixtures. J. Chem. Phys. 1978, 69, 691–696.

(31) Vangara, R.; van Swol, F.; Petsev, D. N. Solvophilic and solvophobic surfaces and non-Coulombic surface interactions in charge regulating electric double layers. J. Chem. Phys. 2018, 148, 044702.

(23) Blum, L.; Henderson, D. Mixtures of hard ions and dipoles against a charged wall: The OrnsteinZernike equation, some exact results, and the mean spherical approximation. J. Chem. Phys. 1981, 74, 1902– 1910.

(32) Lee, J. W.; Nilson, R. H.; Templeton, J. A.; Griffiths, S. K.; Kung, A.; Wong, B. M. 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.

(24) Dong, W.; Rosinberg, M. L.; Perera, A.; Patey, G. N. A theoretical study of the solidelectrolyte solution interface. I. Structure of a hard sphere iondipole mixture near an uncharged hard wall. J. Chem. Phys. 1988, 89, 4994–5009.

(33) Stern, O. Zur theorie der elektrolytischen doppelschicht. Z. Electrochem 1924, 30, 508–516.

(25) Groot, R. D. Density-functional theory for inhomogeneous fluids. Phys. Rev. A 1988, 37, 3456–3464.

(34) Evans, R. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; Marcel Dekker, Inc.: New York, 1992; Chapter 3, pp 85–175.

(26) Perera, L.; Berkowitz, M. L. Dynamics of ion solvation in a Stockmayer fluid. J. Chem. Phys. 1992, 96, 3092–3101.

(35) Wu, J. Density Functional Theory for Chemical Engineering: From Capillarity

ACS Paragon Plus Environment

19

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 24

to Soft Materials. AIChE Journal 2006, 52, 1169–1191.

Spherical Particles. J. Chem. Phys. 2018, 149, 163307.

(36) Rosenfeld, Y. Free-energy model for the inhomogeneous hard-sphere fluid mixture and density functional theory of freezing. Phys. Rev. Lett. 1989, 63, 980–983.

(47) Sandia National Laboratories Tramonto Software: https://software.Sandia.Gov/tramonto/. (48) van Hal, R. E. G.; Eijkel, J. C. T.; Bergveld, P. A general model to describe the electrostatic potential at electrolyte oxide interfaces. Adv. Coll. Interf. Sci. 1996, 69, 31–62.

(37) Roth, R. Fundamental measure theory for hard-sphere mixtures: a review. J. Phys.: Condens. Matter 2010, 22, 063102. (38) Archer, A. J.; Chacko, B.; Evans, R. Spinodal decomposition in a Lennard-Jones fluid. arXiv:1706.08744v1 [cond-mat.soft] 2017,

(49) Tadros, T. F.; Lyklema, J. Adsorption of potential-determining ions at the silicaaqueous electrolyte interface and the role of some cations. J. Electroanal. Chem. 1968, 17, 265–275.

(39) Morse, P.; Feshbach, H. Methods of Theoretical Physics; Feshbach Publishing, 1953; Vol. 1.

(50) Kunz, W.; Henle, J.; Ninham, B. W. The present state of affairs with Hofmeister effects. Curr. Opin. Colloid Interface Sci. 2004, 9, 1–18.

(40) Hill, T. L. An Introduction to Statistical Thermodynamics; Dover: New York, 1986.

(51) W. Kunz, J. H.; Ninham, B. ’Zur Lehre von der Wirkung der Salze’ (about the science of the effect of salts): Franz Hofmeister’s historical papers. Current Opinion in Colloid and Interface Science 2004, 9, 19– 37.

(41) Carnie, S.; Chan, D. Y. C. The structure of electrolytes at charged surfaces: Iondipole mixtures. J. Chem. Phys. 1980, 73, 2949–2957. (42) Ballenegger, V.; Hansen, J.-P. Dielectric permittivity profiles of confined polar fluids. J. Chem. Phys. 2005, 122, 114711.

(52) Lee, S. S.; Fenter, P.; Park, C.; Sturchio, N. C.; Nagy, K. L. Hydrated Cation Speciation at the Muscovite (001)-Water Interface. Langmuir 2010, 26, 16647– 16651.

(43) Gillespie, D.; Nonner, W.; Eisenberg, R. S. Density functional theory of charged, hard-sphere fluids. Phys. Rev. E 2003, 68, 031503.

(53) Graham, D. C. The Electrical Double Layer and the Theory of Electrocapillarity. Chem. Rev. 1947, 41, 441–501.

(44) Barthel, J. M. G.; Krienke, H. Physical Chemistry of Electrolyte Solutions: Modern Aspects; Springer: Darmstadt, 1998.

(54) Frink, L. J. D.; van Swol, F. Solvation forces between rough surfaces. J. Chem. Phys. 1998, 108, 5588–5598.

(45) Valisko, M.; Boda, D.; Gillespie, D. Selective adsorption of ions with different diameter and valence at highly charged interfaces. J. Phys. Chem. C 2007, 111, 15575–15585.

(55) Biesheuvel, P. M.; v. d. Veen, M.; ; Norde, W. A modified Poisson-Boltzmann model including charge regulation for the adsorption of ionizable polyelectrolytes to charged interfaces, applied to lysozyme adsorption on silica. J. Phys. Chem. B 2005, 109, 4172–4180.

(46) Bozic, A. L.; Podgornik, R. Anomalous Multipole Expansion: Charge Regulation of Patchy Inhomogeneously Charged

ACS Paragon Plus Environment

20

Page 21 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(56) L´opez-Garc´ıa, J. J.; Horno, J. PoissonBoltzmann Description of the Electrical Double Layer Including Ion Size Effects. Langmuir 2011, 27, 13970–13974. (57) Plischke, M. Pair correlation functions and density profiles in the primitive model of the electric double layer. J. Chem. Phys. 1988, 88, 2712–2718. (58) Yu, Y. X.; Wu, J. Z.; Gao, G. H. Densityfunctional theory of spherical electric double layers and zeta potentials of colloidal particles in restricted-primitivemodel electrolyte solutions. J. Chem. Phys. 2004, 120, 7223–7233. (59) Tang, Z.; Striven, L. E.; Davis, H. T. Interactions between primitive electrical double layers. J. Chem. Phys. 1992, 97, 9258– 9266. (60) Pizio, O.; Sokolowski, S. On the effects of ion-wall chemical association on the electric double layer: A density functional approach for the restricted primitive model at a charged wall. J. Chem. Phys. 2006, 125, 024512. (61) Jiang, J.; Cao, D.; Henderson, D.; Wu, J. Revisiting density functionals for the primitive model of electric double layers. J. Chem. Phys. 2014, 140, 044714. (62) Blum, L.; Henderson, D. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; Marcel Dekker, Inc.: New York, 1992; Chapter 6, pp 239–276. (63) Grimson, M. J.; Rickayzen, G. Forces between Surfaces in Electrolyte Solutions. Chem. Phys. Lett. 1982, 86, 71–75.

ACS Paragon Plus Environment

21

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Frank van Swol is a Research Professor in the Department of Chemical and Biological Engineering at the University of New Mexico, USA. He received his Ph.D. in physical chemistry from the University of Amsterdam in 1988. He was a Ramsay Memorial Fellow at the Physical Chemistry Laboratory at Oxford University, UK. He then moved to the US, as a visiting faculty at the Chemical Engineering Department of Cornell University and subsequently joined the Chemical Engineering Department at the University of Illinois at Urbana-Champaign. He joined Sandia National Laboratories in 1994 and retired from there in 2015. His research focus concentrates on interfacial phenomena including colloids and electrolytes. He is also interested in additive manufacturing processes such as selective laser melting.

ACS Paragon Plus Environment

22

Page 22 of 24

Page 23 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Dimiter N. Petsev is a professor in the Department of Chemical and Biological Engineering at the University of New Mexico. He received his PhD in physical chemistry at the University of Sofia, Bulgaria in 1996. After graduating from the University of Sofia he spent a year and a half as postdoctoral fellow at Purdue University, followed by research positions in the Center for Microgravity and Materials Research at the University of Alabama in Huntsville, and the Chemical Engineering Department at the University of Houston. His research interests include surface and interfacial phenomena, structure and dynamics of charged colloidal suspensions, soft matter, nanomaterials, electrokinetic phenomena and capillary flows.

ACS Paragon Plus Environment

23

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry

ACS Paragon Plus Environment

24

Page 24 of 24