ARTICLE pubs.acs.org/JPCB
A Classical Density Functional Theory of Ionic Liquids Jan Forsman,*,† Clifford E. Woodward,‡ and Martin Trulsson† † ‡
Theoretical Chemistry, Chemical Centre, P.O. Box 124, S-221 00 Lund, Sweden School of Physical, Environmental and Physical Sciences, University College, University of New South Wales, ADFA Canberra ACT 2600, Australia ABSTRACT: We present a simple, classical density functional approach to the study of simple models of room temperature ionic liquids. Dispersion attractions as well as ion correlation effects and excluded volume packing are taken into account. The oligomeric structure, common to many ionic liquid molecules, is handled by a polymer density functional treatment. The theory is evaluated by comparisons with simulations, with an emphasis on the differential capacitance, an experimentally measurable quantity of significant practical interest.
I. INTRODUCTION Ionic liquids, ILs, have been known to science for almost a century.1 However, they have recently attracted a flurry of activity, mainly due to their potential as alternative, environmentally friendly, solvents. Applications of ILs include specific solvents for heterogeneous and homogeneous catalysis, selective solvents for removal of heavy metal contaminants, electrolytes in supercapacitors, and dispersive agents for nanoparticle stabilization. Some of these applications involve a response to the presence of charged surfaces. For example, many ILs display a characteristic double-humped “camel shape” of the differential capacitance, CD, with a minimum at around zero surface potential. The physical reasons for this behavior are not straightforward and are still being debated in the literature.28,35 Clearly, CD is sensitive to the intermolecular forces and structure of the ILs. For example, in recent work, we showed that attractive dispersion forces in ILs can enhance the depth of the minimum in CD.9 The ability to predict this camel shape has lately become something of a test for a model’s viability for describing ILs . This is especially the case for “coarse-grained” models, of the type we use here, which significantly reduce the number of interaction parameters, compared with so-called all-atom approaches.10 Recent molecular dynamics simulations of CD, using an all-atom model, highlighted the importance of potential parameters and subsequent fluid structure in determining the features of the camel-shaped CD. However, the complexity of the interactions generally makes all-atom models difficult to treat with (semi-) analytic theories, especially for nonuniform fluids. Thus, one is usually constrained to simulations when dealing with these kinds of models. Approximate interaction site integral equation methods have also been used to treat uniform fluids.1113 On the other hand, the accessibilty of coarse-grained models to theories of nonuniform fluids make them particularly attractive. The use r 2011 American Chemical Society
of coarse-grained models also allows the ready identification and prediction of physical mechanisms, which is important to particular phenomena.6 From a theoretician’s viewpoint, ILs are fascinating. In usual aqueous solutions of electrolytes, the concentration is low and the dielectric constant is high. The opposite is true in ILs. These conditions, which lead to crystallization for typical salts, are offset by steric hindrance due to oligomeric nonpolar groups. This has interesting implications on interfacial phenomena. In typical electrolytes, the structures at charged interfaces are broadly determined by the Debye length. In ILs, the steric length-scale is also relevant. Thus, ILs provide us with a unique opportunity to investigate a system with strong interplay between electrostatic and steric interactions. Such systems cannot be accurately described by a simple mean-field theory, such as the PoissonBoltzmann approximation. Instead, one requires a theory that is able to describe electrostatic and steric correlations. Density functional theory, DFT, provides a flexible and versatile framework within which systematic improvements to the PoissonBoltzmann approximation can be incorporated to better describe these mechanisms. In particular, one is able to include ionion correlations and steric effects, so important in this class of fluids. As pointed out by Forsman,14 the length scales for steric and electrostatic correlations may be quite different. However, as stated above, we anticipate that in ILs, the electrostatic and steric length scales are comparable. In this work, we borrow ideas from previous studies on the one-component plasma (OCP)1517 and simple salt solutions14 to construct a DFT wherein charge correlations are Received: December 10, 2010 Revised: January 28, 2011 Published: April 01, 2011 4606
dx.doi.org/10.1021/jp111747w | J. Phys. Chem. B 2011, 115, 4606–4612
The Journal of Physical Chemistry B accounted for. The resultant theory is phenomenological in nature but, as we will see, agrees quite well with simulations. Electrostatic correlations have been shown to be important in interfacial phenomena of highly coupled electrostatic fluids. As with dispersion forces, electrostatic correlations provide an attractive component to the fluid energy. Thus, we expect that they will also play a role in establishing the camel shape of the CD. That this is the case is seen from both simulations18,19 and modified Poisson Boltzmann treatments20 of a restricted primitive model (RPM) electrolyte. These investigations show that the RPM displays a camel-shaped CD at low temperatures and densities. However, electrostatic correlations do not tell the whole story for ILs. That is, the CD becomes bell-shaped (i.e., possesses a single maximum) for the RPM at densities and temperatures typical for ILs. This suggests that the simple RPM is an inadequate model for ILs and other factors, such as dispersion forces, and molecular structure8,9 plays an important role in determining their properties. In this work, we will use a model for ILs wherein the molecular ions interact via dispersion, as well as electrostatic, forces and the cation has an oligomeric structure, consistent with imidazolium ions. Thus, this model also requires us to use the methods of polymer DFT.21 The result will be a simple, albeit approximate, density functional theory for ILs, IL-DFT. In some respects, the IL-DFT can be seen as a logical progression from recent approximate theories for ILs developed by Kornyshev et al.6 and Lauw et al.7,22 In the former case, steric effects were approximated by use of a lattice model. The latter model, also based on a lattice model, included oligomeric structures and a localized dielectric response, related to the polarizability of separate parts of the molecules. Both of these theories were able to demonstrate the existence of a camel-shaped CD, under some conditions. The IL-DFT is an off-lattice theory whose accuracy can be straightforwardly tested against computer simulations of the same molecular model. Here, we will test the IL-DFT against our earlier Monte Carlo simulations for CD.9 The prediction of a camel-shaped CD at appropriate values for density and temperature is good support for a theory purporting to describe, at least qualitatively, the behavior of ILs.
II. MODEL Recently, we presented a computer simulation study of ILs at charged interfaces.9 We shall use the same model in this work. The cations are modeled as pentamers composed of bonded Lennard-Jones (L-J) monomers. The Lennard-Jones pair potential, ΦL-J, is given by " 6 # σ 12 σ ΦL-J ðrÞ ¼ 4εL-J ð1Þ r r where r is the separation between particles. The size parameter, σ, also defines the distance between bonded monomers. The bond potential, Vb, contains no angular part and is most easily represented by the Boltzmann factor: eβVb(rij) µ δ(rij σ), where rij is the distance between connected monomers i and j and β = 1/kT is the inverse thermal energy. Thus, the pentamers are fully flexible, forming a “pearl-necklace” of L-J particles. The second monomer in the chain carries a unit of positive charge. The anions are modeled as simple monomeric particles, which interact via the same L-J parameters as the monomers of the cation. All simulations were performed with εL-J = 0.6kT. Two choices of σ were considered, 4 and 6 Å, and temperature was set at 393
ARTICLE
K.4 The electrostatic interaction between two charges, qR and qλ, separated by r is given by ΦRβ el ðrÞ ¼
1 qR qλ 4πε0 εr r
ð2Þ
where ε0 is the permittivity of vacuum, and εr is due to the contributions to the low frequency dielectric response additional to those not already within our model. Typical experimental values for the low-frequency dielectric constant of ILs range between 10 and 15. Because our model is a “coarse-grained” one, dielectric effects due to electronic, orientational, and distortional polarization are unaccounted for. The electronic contribution to εr is around 2. The latter two contributions are due to charge distribution and bond length and angle flexibility, respectively. In our model, the molecules possess no permanent dipoles with respect to the center of charge and, compared with “real” ILs, the electrostatic coupling may be overestimated. Recent work by Izgordina et al.23 show that the orientational contribution to εr, can be as large as 20, depending on the size of the molecular dipole moments, while the distortion contribution remains small (approximately unity). The density functional theory that we will use in this work is based on an OCP model, wherein the important ion correlation effects are assumed to be short-ranged, that is, comparable to the range of steric effects. Hence, the long-range dielectric screening of ions will also not be accounted for. This is the so-called ionic polarization contribution to the dielectric constant.23 Thus, we have used a typical experimental value, εr = 12, in the DFT calculations. The response of the differential capacitance to changes of εr was investigated in our previous work.9 We recall that the purpose of the current work is not to fine-tune parameters to specifically describe a particular IL. On the other hand, simulations do capture a low-frequency ionic polarization contribution to εr. The appropriate value for the dielectric constant to be used in simple simulation models such as ours is subject to considerable debate in the literature.2325 Here, we have used ε*r = 7.6 in all our simulations. Note that we have used an asterix in our nomenclature to differentiate the simulation dielectric constant. Thus, in our simulations, Coulomb interactions are scaled down by 1/7.6 relative to the free space value. The electrodes are modeled as uniformly charged, infinite planar surfaces normal to the z axis. One surface has a charge density σel; the other carries a charge density of equal magnitude but opposite sign. All fluid particles also experience an underlying, nonelectrostatic interaction with the electrodes, obtained by integrating the Lennard-Jones potential over the half-spaces constituting the surfaces, " # 2 σw 9 1 σw 3 aw ð3Þ βwL-J ðzÞ ¼ 2πεL-J 45 z 3 z where z is the distance from the plane of charge. In most cases, we have set σw = σ, but we also investigated the effect of varying σw/ σ. The adsorption parameter aw is either 0 or 1, modeling either an underlying nonadsorbing or adsorbing surface interaction, respectively. The total nonelectrostatic interaction due to the surfaces, VL-J, is given by the sum of the two separate contributions, that is, VL-J = wL-J(z) þ wL-J(h z), where h is the distance between the surfaces. Simulation Method. Monte Carlo simulations were performed in the Canonical Ensemble. We chose a surface separation large 4607
dx.doi.org/10.1021/jp111747w |J. Phys. Chem. B 2011, 115, 4606–4612
The Journal of Physical Chemistry B
ARTICLE
Figure 1. Monomer density distributions from canonical simulations, with neutral nonadsorbing (aw = 0) surfaces. The number of salt pairs, Ns, is different for each curve, as indicated. T is the temperature; Fm is the total monomer density. (a) T = 393 K. Here, we observe a response typical to a liquid —that is, a liquidgas interface is formed—to accommodate the reduced average density. (b) T = 786 K. A typical response by a supercritical fluid, where the density drops uniformly, throughout the slit.
enough to generate bulklike conditions near the midplane of the resulting slit. The simulations required some preliminary testing to establish the prescribed bulk density at the central region of the simulated volume. Periodic boundary conditions were implemented in the xy plane, parallel with the surfaces. A mean-field correction, calculated using the average density profiles within the simulation volume, accounted for electrostatic interactions originating from outside the simulation box. The shorter-ranged dispersion terms were treated using the minimum image convention. In previous work,9 we demonstrated that this kind of treatment of boundary conditions agrees well with grand canonical ensemble simulations, combined with Ewald sums (GCMC þ Ewald) for the long-ranged corrections. For one of the cases in this study (aw = 1), we also performed a comparison with the GCMC þ Ewald approach. In this work, we wish to focus on the differential capacitance, CD, because it provides a sensitive measure of the interplay between electrostatic and steric mechanisms. It is given by the formula: CD = dσel/dΨS, where ΨS is the average potential at the surface, with charge density σel. The average potential at any point z, Ψ(z), was determined relative to the midplane value, where bulk conditions were assumed, that is, Ψ(h/2) = 0. The differential capacitance was obtained as a numerical derivative. We used the following approach to determine whether our model describes a proper liquid, that is, that we are below the critical temperature under the conditions of our simulation. Canonical Ensemble simulations were performed for the fluid in the presence of neutral, nonadsorbing surfaces (aw = 0, σs = 0). The number of salt pairs in the fluid were gradually decreased during the simulation. For a liquid phase, this process will lead to dewetting at one or both of the surfaces and the buildup of at least one gasliquid boundary. On the other hand, above the critical temperature, the process will generate a gradual lowering of the entire density profile. Both of these behaviors are illustrated in Figure 1.
III. IL-DFT It is convenient to develop the theory in three stages: (i) We first consider a “polymer solution” of L-J pentamers and monomers in the presence of the L-J surfaces as described above. (ii) Second, the charges of the ions and the surfaces are switched on, and their effect is modeled using a meanfield (PB) approximation. (iii) Contributions to the free energy due to charge correlations are added.
The L-J polymer solution is treated using the polymer-DFT.21 Because this theory has become “standard” in recent years, we shall only provide an outline here. A more detailed account of the theory can be found elsewhere.21,26 We consider a fluid of linear oligomers, each with r monomers, labeled by the subscript “m”. In the IL model considered here, these oligomers have r = 5 and will become cations upon switching on the charges. In addition, there are “solvent” particles, with index “s”. These will become the anions. The solvent particles and monomers have the same L-J diameter, σ, although the theory has been developed to treat differently sized particles, as well.27,26 We shall invoke the commonly used approximation that the repulsive part of the L-J is modeled as a hard core, which acts at a particleparticle separation of σ. The attractive part is handled in a mean-field fashion. A chain configuration is represented as R = (r1, ..., rr), where ri is the coordinate of monomer i. The Q total bonding potential is given by, VB(R), i.e. eβVB(R) µ r1 i=1 δ(|riþ1ri|σ). It is convenient to define a density distribution N(R), such that N(R) dR is the number of polymer molecules having configurations between R and R þ dR. For a fluid of ideal chains, the exact 21 expression for the free energy functional, F id p , is Z id βF p ½NðRÞ ¼ NðRÞðln½NðRÞ 1Þ dR þ Z β
NðRÞ VB ðRÞ þ Vex ðRÞ dR
ð4Þ
where Vex(R) is an external potential, in this case due to the surfaces. The total free energy will contain ideal terms from the solvent as well as excess contributions due to the steric hard sphere and attractive L-J interactions between particles. In the presence of L-J surfaces, a general expression for the grand potential, ΩL-J, of the L-J polymer fluid, can be written as βΩL-J NðRÞ, Fs ðrÞ ¼ βF idp ½NðRÞ þ
Z
Fs ðrÞ ln Fs ðrÞ 1 dr
Z þ βF hs Fm ðrÞ, Fs ðrÞ þ β VL-J ðrÞ μp Fm ðrÞ dr Z þ β VL-J ðrÞ μs Fs ðrÞ dr þ 4608
β 2
Z Z
0
∑m ∑s Fm ðrÞ Fs ðr0Þ ΦL-J ðjr r0jÞ dr dr0
ð5Þ
dx.doi.org/10.1021/jp111747w |J. Phys. Chem. B 2011, 115, 4606–4612
The Journal of Physical Chemistry B
ARTICLE
where μp and μs are the polymer and “solvent” chemical potentials, respectively. The prime above the second integral for the L-J interactions, indicates that regions where |r r0 | e σ should be excluded. The hard sphere, exclusion-free energy, Fhs, 26,27 arrived at has to be approximated. Forsman and co-workers an expression by integrating an equation of state developed by Wichert et al.28 The expressions for the weighted29 densities, Fm(r) and Fs(r), are in our case very simple because they have identical diameters: Z 3 Fm ðrÞ ¼ F ðr0 Þ dr0 ð6Þ 4πσ3 jr r0 j < σ m Upon switching on the electrostatic interactions, an additional term must be added to ΩL-J so as to obtain the total grand potential Ωtot. We let Ωel denote this extra term. To provide an intermediate reference, Ωel is initially approximated at the meanfield (i.e. PoissonBoltzmann) level. We denote this reference functional by IL-PB, with Z Z β 0 0 nR ðrÞ nλ ðr0 Þ ΦRλ βΩel ½IL-PB ¼ el ðjr r jÞ dr dr 2 R λ Z þβ Vel ðrÞ nR ðrÞ dr ð7Þ
∑∑
∑R
Here, nR(r) are ionic densities, where the indexes R and λ are summed over all the charged species. Note that for anions, the ionic and solvent densities, Fs(r) (see above), are the same quantities, but for cations, the ionic density corresponds to the density distribution of the second monomer in the chain, that is, the charged monomer. The first term in eq 7 is just the mean-field electrostatic interaction between ions. The second term describes the interactions between the ions and the uniform surface charge. We note that Ωel[IL-PB] is zero for the bulk phase IL. The total free energy functional obtained as Ωtot[IL-PB] = ΩL-J þ Ωel[IL-PB], and represents an intermediate, mean-field version of our theory, which corrects the traditional PB approximation by accounting for dispersion forces. This level of theory is an “off-lattice” analogue of the steric lattice model recently described by Kornyshev,6 but generalized to include attractive dispersion forces. As shown by Kornyshev, the adsorption saturation induced by steric effects causes a decrease in CD at high surfaces charges. Furthermore, in earlier work,9 we showed that the attractive dispersion interactions enhance the minimum in CD at low surface charge. Thus, we expect that the IL-PB approximation will predict a camel-shaped CD. This notwithstanding, the IL-PB level of theory ignores electrostatic correlations, which provide another source for lowering the free energy of the fluid. A significant result in the current work is the introduction of an additional free energy term to Ωel[IL-PB] to account for electrostatic correlations. The resulting functional will be denoted IL-DFT and has the following form: Z
βΩel ½IL-DFT ¼ βΩel ½IL-PB þ β
∑R nR ðrÞ
fc ½^nR ðrÞ μRc dr ð8Þ
where μRc is obtained by differentiating fc at the bulk density of R. Electrostatic correlations, together with dispersion interactions, provide a cohesive binding energy to the ionic liquid, favoring the liquid phase. This additional binding energy is tempered by the loss in entropy that always accompanies correlations. The strength of the electrostatic forces in ILs will mean that
correlations are an important contribution to the free energy. As explicitly demonstrated in previous work,14,30 the most significant effect of these correlations occurs because of the mutual repulsion between like-charged species. These can be incorporated into the DFT using an idea recently suggested by Forsman14 wherein correlations between ions with the same charge can be approximated using the OCP model. Here, we shall use the “hole corrected DebyeH€uckel” theory (DHH) developed by Nordholm and co-workers1517 for the OCP. This theory provides us with a natural length scale over which to estimate the correlational free energy in a manner analogous to that used for steric effects in the L-J fluid described above. In the bulk OCP, the mean-field energy is zero. On the other hand, correlation effects are due to the mutual repulsion between particles of the plasma in the presence of a uniform neutralizing background. The electrostatic repulsion means that particles in the OCP exclude each other so as to create an electrostatic correlational hole of radius hc. This hole has a minimum value equal to the steric radius, σ, and is generally larger than that in dilute electrolytes. In ILs, however, we expect hc ≈ σ. The distance hc provides a natural length for correlational effects in the nonuniform IL. Thus, we define the “coarse-grained” ionic density, n^R: Z 3 ^nR ðrÞ ¼ nR ðr0 Þ dr0 ð9Þ 4πhc 3 jr r0 j < hc In the DHH treatment, a closed analytic expression can be obtained for hc and the corresponding free energy. For a bulk density equal to n^R(r), the free energy per particle of the OCP is16 1 2π βfc ½^nR ðrÞ ¼ 1 þ pffiffiffi ln½3 þ ln½ω2 þ ω þ 1 ω2 4 3 3 2 2ω þ 1 pffiffiffi tan1 pffiffiffi ð10Þ 3 3 which is used in the last term of eq 8. Here, h i1=3 ω ¼ 1 þ c^nR 1=2 with
3 rffiffiffiffiffiffi
q 4π
R 3=2 c ¼ ð3lB Þ
e 3
ð11Þ
ð12Þ
Here, lB is the Bjerrum length, lB = e2β/4πε0εr. In this work, we shall make the simple approximation that hc is determined by the ion density in the bulk, nbR,16 hc ¼
ω1 k
ð13Þ
where κ = (4πlB^nbR)1/2 is the Debye screening length in the bulk. Note that at high densities, there are some corrections due to steric effects,17 which have been neglected here. This should be a reasonable approximation in our case, given that hc ≈ σ for the systems we have investigated. As discussed earlier, the correlation correction used here accounts for only short-ranged effects, expected to dominate in ILs, and does not account for the ionic polarization contribution to the dielectric screening of electrostatic interactions. Thus, the mean-field terms (i.e. the first two terms in eq 8) remain unaffected by the correlational term, and 4609
dx.doi.org/10.1021/jp111747w |J. Phys. Chem. B 2011, 115, 4606–4612
The Journal of Physical Chemistry B
ARTICLE
Figure 2. Bulk pressure vs the sum of monomer and simple ion densities, for our two different IL-DFT approaches, with σ = 6 Å, are given by thick, black and blue lines. In IL-DFT(ε) calculations, εL-J is equal to the value adopted in the simulations, that is, βεL-J = 0.6. This parameter is established in a different manner for the IL-DFT(F) version, as explained in the text. In this particular case, this leads to βεL-J = 0.91. The corresponding curve for our reference IL-PB functional, which lacks ion correlation contributions, is shown by a thin, violet line. The difference between the IL-PB(ε) and the IL-DFT(ε) curves is entirely due to ion correlations in the bulk fluid. The densities at which calculations of the differential capacitance were made are indicated by short, horizontal lines. The arrow indicates the (liquid) density at which there is gasliquid coexistence in the simulations (cf. Figure 1).
the ionic polarization contribution must be included in the value of εr. The total grand potential, Ωtot[IL-DFT] is obtained as Ωtot ½IL-DFT ¼ ΩL-J þ Ωel ½IL-DFT
ð14Þ
The functionals described above can be simplified by making use of the flat surface geometry, whereby the (x, y) dimensions, parallel with the “electrode walls”, can be integrated out. Minimization of the chosen functional, with respect to the fluid particle densities then follows in the usual manner, wherein electroneutrality is ensured via a Donnan potential. In the results presented below, our focus will, of course, be on the IL-DFT level of theory. However, we will occasionally include IL-PB results to highlight effects from ion correlations. A. Mapping for Simulation Comparisons. The approximate nature of the IL-DFT introduces a degree of uncertainty as to how comparisons with the simulations should be attempted. As noted earlier, we set the dielectric constant, εr, to 12 in the ILDFT calculations, a value similar to those typically found in experiments. We also chose the parameter σ to be the same as in the simulations. On the other hand, varying the choice of the L-J strength parameter, εL-J, will give rise to different values for the total bulk density of monomers and simple ions, Fm(bulk) þ Fs(bulk). IL-DFT and simulations will generally not agree on the value of this bulk density for a given choice of εL-J. A consequence of this is that the IL-DFT bulk phase diagrams will will be quantitatively different from the corresponding simulated ones. From our earlier discussions, it is clear that the simulation is performed at an average liquid density close to, but slightly greater than, that at liquidgas coexistence; that is, we are simulating a stable liquid phase. We need some procedure with which to establish the “corresponding” fluid states for comparison with the IL-DFT calculations. We considered two different choices:
Figure 3. Simulated and calculated differential capacitance curves, with aw = 0 (nonadsorbing) and σ = 6 Å.
• Simply set the εL-J used in the IL-DFT calculations be equal to the value adopted in the simulations; that is, βεL-J = 0.6. We shall denote this approach IL-DFT(ε). • Choose εL-J in the IL-DFT such that the same liquid coexistence density is obtained, as in the simulations. We denote this version IL-DFT(F). The advantage of the first choice, IL-DFT(ε), is that it is easily implemented. A disadvantage is that, in principle, it is possible that the simulations one wants to compare with are performed at a temperature that is above the critical point predicted by the ILDFT. With the inclusion of the correlation term, eq 10, we never encountered the latter problem; that is, the IL-DFT bulk equation of state is then sufficiently accurate for the parameter set investigated. However, without this term (i.e., with the ILPB), the system just barely displays bulk demixing at βεL-J = 0.6, as illustrated below. To ensure a stable liquid phase in the ILDFT, we used a pressure 10 times higher than at the gas phase spinodal, where the pressure versus density has a local maximum. This approach was adopted for all DFT calculations, including IL-PB. Because the liquid has a low compressibility, we note that the corresponding density is only slightly higher than that at coexistence. The resulting pressure-versus-density curves for the IL-DFT using the two different approaches, (ε) and (F), are illustrated in Figure 2. We see that the bulk densities become quite different, depending on which simulation comparison approach we use. However, as demonstrated below, this has a remarkably small impact on the predicted differential capacitance curve. Also included in Figure 2 is the corresponding IL-PB(ε), which results from removing the correlation contribution to IL-DFT(ε). In this case, the system is very close to the critical point, and we clearly see that the correlations make a substantial contribution.
IV. RESULTS The accuracy of the IL-DFT will be evaluated via comparisons between simulated and predicted values for the differential capacitance, CD(ΨS), as a function of surface potential. In Figure 3, such curves are presented for a model with an underlying repulsive surface, aw = 0, and σ = 6 Å. The agreement is qualitatively excellent, if not quantitative. For inherently repulsive surfaces, depletion occurs near the surfaces at low surface potentials. This depletion comes from several sources. Dispersion and electrostatic correlations lead to a positive surface energy contribution due to both ionic species, leading to a diminished concentration at the surfaces. Furthermore, the 4610
dx.doi.org/10.1021/jp111747w |J. Phys. Chem. B 2011, 115, 4606–4612
The Journal of Physical Chemistry B
Figure 4. Simulated and calculated differential capacitance curves, with aw = 0 (nonadsorbing) and σ = 4 Å. IL-DFT(F) data are not shown, but they agree closely with the displayed IL-DFT(ε) curve. We have also included predictions by the reference IL-PB functional, which neglects ion correlations.
Figure 5. Simulated and calculated net charge density profiles at a negatively charged nonadsorbing electrode, with σ = 4 Å and 1/σs = 340 Å2/e. The net charge density is scaled by the bulk salt concentration, nbþ. Results from IL-DFT(ε) calculations are shown, but also from our simpler reference functional, IL-PB(ε), that lacks contributions from ion correlations.
restriction of the configurational entropy of the cationic oligomers also causes depletion from repulsive surfaces. Small changes in the surface charge density require relatively large changes in the surface potential, and the corresponding differential capacitance is low. At intermediate absolute values of the surface potential, the screening of the charges increases, and the response of the surface potential to changes in the surface charge density, σel, starts to saturate, leading to an increasing CD. At very high surface potentials, on the other hand, steric effects start to substantially counteract an increased density of the adsorbed ions. Instead, they adsorb in layers, which leads to large changes in the surface potential with respect to σel and a consequent decrease of CD. These mechanisms account for the so-called “camel hump” shape that is commonly observed with ILs.28 The simulated CD(ΨS) curve in Figure 3 has this characteristic shape, which is also reasonably accurately reproduced by the IL-DFT, although steric effects seem to set in at higher values of ΨS than observed in simulations. In Figure 4, we show a comparison of CD calculated using ILDFT(ε) and simulations, but with σ = 4 Å. Here, we see a deeper and wider minimum in the simulated CD, which is qualitatively reproduced by the IL-DFT. To reduce the number of curves, we have not included results by the IL-DFT(F), but the agreement with IL-DFT(ε) is even closer than was observed in Figure 3. The
ARTICLE
Figure 6. Simulated and calculated differential capacitance curves, with aw = 1 (adsorbing) and σ = 6 Å.
simulated values appear somewhat more “U-shaped” than the density functional predictions, but the agreement seems excellent, considering the simplicity of the theory. Although not seen on this scale, CD(Ψs) will reach a maximum at a high positive surface potential. Specifically, this maximum is located at ∼1.5 V. Also shown in Figure 4 are CD predictions from the IL-PB functional. Although the same qualitative features are observed, the quantitative performance is severely deteriorated, as compared with the IL-DFT functional. Clearly, the free energy lowering due to electrostatic correlations and the interaction with the surfaces has a significant effect on CD, which is largely captured by our simple functional. To explore the effect of correlations further, we compare the charge density profiles obtained at the different levels of theory. In Figure 5, we display net charge density profiles outside a negatively charged surface of intermediate strength, as obtained via simulations, and IL-DFT(ε) as well as IL-PB(ε) calculations. The IL-DFT(ε) predictions agree qualitatively with simulation data, reproducing the oscillatory behavior as well as a small “hump” at z = 78 Å. On the other hand, the IL-PB(ε) results significantly overestimate the counter-charge peak adjacent to the surface and completely fails to capture the oscillatory behavior found by simulations. This early onset of charge saturation is the reason behind the underestimation of the minimum in CD by the IL-PB functional. Though closer to the simulations, the quantitative performance of the IL-DFT level of theory is still far from perfect, with a too broad primary peak in the charge density, which results in a “phase lag” of the oscillations further out. Still, it should be noted that this initial and rather crude estimate of ion correlations does lead to substantial improvements relative to the IL-PB. There are several options for improvement of the IL-DFT. These include the introduction of other weight function alternatives for the coarse-grained densities,31,32 a variationally determined “Coulomb hole”,33 and an improved estimate of how excluded volume influences the interactions. At least some of these options will be explored in future work, where we also will consider the introduction of chain stiffness34 into the model. We have noted that the minimum of CD(ΨS) develops as a result of density depletion at near-neutral nonadsorbing electrodes. But what if these electrodes adsorb the IL molecules, even when there are no surface charges? In that case, one would anticipate that the camel hump is reduced or even disappears. This conjecture was, indeed, corroborated in a previous simulation study.9 Similar results are found in Figure 6, together with corresponding IL-DFT predictions. Again, we find satisfactory 4611
dx.doi.org/10.1021/jp111747w |J. Phys. Chem. B 2011, 115, 4606–4612
The Journal of Physical Chemistry B agreement, which seems to capture all essential features of this system with a high accuracy. The simulation data are a bit more noisy with adsorbing electrodes, particularly at low surface potentials. We note that the IL-DFT curves display a tendency to form a minimum, but that the nonelectrostatic adsorption is strong enough to prevent this from actually forming.
V. SUMMARY We have presented a simple and versatile yet remarkably accurate IL-DFT . This was achieved by combining standard polymer-DFT with a correlation-corrected DFT for salt solutions. In the latter case, ideas from the hole-corrected Debye H€uckel theory of one-component plasmas are borrowed to account for ionic correlations. The inclusion of these correlations generates a theory that takes into account that there is a net electrostatic cohesion in the bulk liquid, which is not recognized at the level of the PoissonBoltzmann approximation. The present theory lends itself to further improvements. For instance, polymer-DFT has been developed to take intrinsic chain stiffness into account34 or differently sized particles.26,27 These developments can be directly imported into the present theory. ’ AUTHOR INFORMATION
ARTICLE
(19) Vatamanu, J.; Borodin, O.; Smith, G. D. J. Am. Chem. Soc. 2010, 132, 14825. (20) Lamperski, S.; Outhwaite, C.; Bhuiyan, L. J. Phys. Chem. 2009, 113, 8925. (21) Woodward, C. E. J. Chem. Phys. 1991, 94, 3183. (22) Lauw, Y.; Horne, M. D.; Rodopoulos, T.; Nelson, A.; Leermakers, F. A. M. J. Phys. Chem. B 2010, 114, 11149. (23) Izgordina, E. I.; Forsyth, D. R.; MacFarlane, M. Phys. Chem. Chem. Phys. 2009, 11, 2452. (24) Krossing, I.; Slattery, C.; Daguenet, J. M.; Dyson, P. J.; Oleinikova, A.; Weing€artner, H. J. Phys. Chem. B 2006, 110, 12682. (25) Kobrak, M. N.; Li, H. Phys. Chem. Chem. Phys. 2010, 12, 1922. (26) Forsman, J.; Woodward, C. E. J. Chem. Phys. 2004, 120, 506. (27) Forsman, J.; Woodward, C. E.; Freasier, B. C. J. Chem. Phys. 2002, 117, 1915. (28) Wichert, J. M.; Gulati, H. S.; Hall, C. K. J. Chem. Phys. 1996, 105, 7669. (29) Nordholm, S.; Johnson, M.; Freasier, B. C. Aust. J. Chem. 1980, 33, 2139. (30) Forsman, J. J. Phys. Chem. B 2004, 108, 9236. (31) Tarazona, P. Phys Rev. A. 1985, 32, 2672. (32) Rosenfeld, Y. Phys. Rev. Lett. 1989, 63, 980. (33) Hooper, M. A.; Nordholm, S. Aust. J. Chem. 1981, 34, 1809. (34) Forsman, J.; Woodward, C. E. J. Chem. Phys. 2003, 119, 1889. (35) Fedorov, M. V.; Kornyshev, A. A. J. Phys. Chem. B 2008, 112, 11868.
Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT J.F. acknowledges financial support by the Swedish Research Council. ’ REFERENCES (1) Walden, P. Bull. Acad. Sci. St. Petersburg 1914, 1, 405. (2) Nanjundiah, C.; McDevitt, S. F.; Koch, V. R. J. Electrochem. Soc. 1997, 144, 3392. (3) Alam, M. T.; Islam, M. M.; Okajima, T.; Oshaka, T. Electrochem. Commun. 2007, 9, 2370. (4) Lockett, V.; Sedev, R.; Ralston, J.; Horne, M.; Rodopoulos, T. J. Phys. Chem. C 2008, 112, 7486. (5) Alam, M. T.; Islam, M. M.; Okajima, T.; Oshaka, T. J. Phys. Chem. C 2008, 112, 16568. (6) Kornyshev, A. A. J. Phys. Chem. B 2007, 111, 5545. (7) Lauw, Y.; Horne, M. D.; Rodopoulos, T.; Leermakers, F. A. M. J. Phys. Chem. B 2009, 111, 5545. (8) Fedorov, M. V.; Georgi, N.; Kornyshev, A. A. Electrochem. Commun. 2010, 12, 296. (9) Trulsson, M.; Algotsson, J.; Forsman, J.; Woodward, C. E. J. Phys. Chem. Lett. 2010, 1, 1191. (10) Lynden-Bell, M. R.; Del Papolo, M. G.; Tristan, G. A.; Kohanoff, J.; Hanke, C. G.; Harper, J. B.; Pinilla, C. C. Acc. Chem. Res. 2007, 40, 1138. (11) Bruzzone, S.; Malvaldi, M.; Chiappe, C. Phys. Chem. Chem. Phys. 2007, 9, 5576. (12) Bruzzone, S.; Malvaldi, M.; Chiappe, C. J. Chem. Phys. 2009, 129, 074509. (13) Malvaldi, M.; Bruzzone, S.; Chiappe, C.; Gusarov, G.; Kovalenko, A. J. Phys. Chem. B 2009, 113, 3536. (14) Forsman, J. J. Chem. Phys. 2009, 130, 064901. (15) Nordholm, S. Chem. Phys. Lett. 1984, 105, 302. (16) Penfold, R.; Nordholm, S.; J€onsson, B.; Woodward, C. E. J. Chem. Phys. 1990, 92, 1915. (17) Penfold, R.; Nordholm, S. J. Chem. Phys. 1991, 95, 2048. (18) Lamperski, S.; Zydor, A. Electrochim. Acta 2007, 52, 2429. 4612
dx.doi.org/10.1021/jp111747w |J. Phys. Chem. B 2011, 115, 4606–4612
