11868
2008, 112, 11868–11872 Published on Web 08/26/2008
Ionic Liquid Near a Charged Wall: Structure and Capacitance of Electrical Double Layer Maxim V. Fedorov*,†,‡ and Alexei A. Kornyshev*,§ Max Planck Institute for Mathematics in the Sciences, D 04103, Leipzig, Germany, UnileVer Centre for Molecular Science Informatics, Department of Chemistry, UniVersity of Cambridge, Lensfield Road, Cambridge, CB2 1EW, U.K., and Department of Chemistry, Faculty of Natural Sciences, Imperial College London, SW7 2AZ London ReceiVed: April 21, 2008; ReVised Manuscript ReceiVed: July 12, 2008
We study the effects of ion size asymmetry and short-range correlations on the electrical double layer in ionic liquids: we perform molecular dynamics simulations of a model ionic liquid between two “electrodes” and calculate the differential capacitance of each as a function of the electrode potential. The capacitance curve has an asymmetric “bell-shape” character, in qualitative agreement with recent experiments and the mean-field theory (MFT) which takes into account the limitation on the maximal local density of ions. The short-range ionic correlations, not included in the MFT, lead to an oVerscreening effect which changes radically the structure of the double layer at small and moderate charging. With the radius of cations taken to be twice as large as anions, the position of the main capacitance maximum is shifted positively from the potential of zero charge (PZC), as predicted by MFT. An extension of the theory (EMFT), however, reproduces the simulated capacitance curve almost quantitatively. Capacitance curves for real ionic liquids will be affected by nonspherical shape of ions and sophisticated pair potentials, varying from liquid to liquid. But understanding the capacitance behavior of such model system is a basis for rationalizing those more specific features. 1. Introduction The wave of interest to room temperature ionic liquids (ILs)1 has recently reached electrochemistry and its applications.2 A feature article in this journal, focused on fundamentals of the double layer theory in ILs,3 was followed by publications reporting new capacitance measurements for different electrodes in a wide range of electrode potentials.4 Reference 3 stressed that one cannot apply the Gouy-Chapman-Stern theory to these dense ionic systems and discussed other approaches. As the simplest one, an alternative meanfield theory (MFT) was suggested that took into account constraints on the ion packing in ILs, the so-called lattice saturation effect. That approach resulted in a formula for the diffuse double layer capacitance as a function of electrode potential. (Similar expressions were independently obtained by several authors—mainly in the context of concentrated electrolytic solutions.5) The Gouy-Chapman formula followed from it as a “limiting case” but never realizable for dense ILs. Furthermore, ref 3 extended the MFT to describe ILs with different sizes of cations and anions. Here we report the results of molecular dynamics (MD) simulations for the capacitance response of simple dense IL made of positively and negatively charged Lennard-Jones spheres. Studying this model system is a necessary step before rationalizing the double layer properties of more complex ILs. * Corresponding authors. E-mail:
[email protected],
[email protected] (M.V.F.);
[email protected] (A.A.K.). † Max Planck Institute for Mathematics in the Sciences. ‡ University of Cambridge. § Imperial College London.
10.1021/jp803440q CCC: $40.75
This study is targeted to test the predictions of the mean-field theory and answer the following questions:3 1. How does the double layer capacitance curve look? 2. What are the counterparts of the obtained response? 3. What is the role of the asymmetry of the ion sizes? 4. How good does the extended MFT (EMFT) describe the results? 2. Model and Methods The IL was modeled as a 1 to 1 mixture of counterlike singly charged spheres with a short-range repulsive Lennard-Jones ij potential uLJ (r) ) 2kBT[(ri0 + rj0)/r]12 between the particles i and j (kB ) Boltzmann constant, T ) temperature). The radii of the spheres r0 were 0.5 nm for cations and 0.25 nm for anions, reflecting an asymmetry of ion sizes in ILs.6 We screened all Coulomb interactions in the system by an effective dielectric constant ε* ) 2.0, which accounts for electronic polarizability of the particles. We put 1050 cations and 1050 anions between two electrodes in a periodic rectangular box with lengths in the X and Y directions equal to 11 nm and the length in the Z direction equal to 40 nm. The electrodes were modeled as two parallel XY square lattices of densely packed Lennard-Jones spheres with radii 0.11 nm; surface charge densities on the electrodes were varied by the partial charge of those spheres. Each pair of electrodes was separated by the 24 nm distance filled by IL, with a 16 nm slab of vacuum separating it from the next periodic image in Z direction. For simulations we used Gromacs 3.3 software.7 The electrostatic interactions were treated with use of particle-mesh Ewald summation, corrected for slab geometry.8 We performed 2008 American Chemical Society
Letters
J. Phys. Chem. B, Vol. 112, No. 38, 2008 11869
Figure 1. Differential capacitance, C, as a function of electrode potential U. Points shown as squares are calculated from the results of molecular dynamic simulations (as described in Models and Methods and in the Supporting Information) performed for a model IL in which cations are twice as large as anions; to guide eye, they are connected by dashed lines. A solid line displays the results of EMFT calculation of the capacitance of such system (for details and parameters, see ref 14). The potential is calculated relative to the bulk; the potential of zero charge in this model corresponds to U ≈ 0.14 Thus, the main maximum of capacitance is shifted to the right of PZC. The arrow depicts the position of the secondary maximum on the simulation curve in the region of negative potentials. The comparison of the two curves reveals the expected agreement at the wings, and some “fine structure” features at moderate potentials that are missing in EMFT, but the basic qualitative agreement between simulations and the EMFT is astounding. The inset on the right displays (with the permission of the authors of ref 17) the experimental data on the differential capacitance of hmimCl at 372 K.17 Detailed comparison of the simulations with these data does not make sense due to the simplicity of the simulation model, but the main features in the experimental curve remind us those seen in the simulation curve. The inset on the left answers the question how well the asymptotic laws at the wings are reproduced by the simulations, showing the results in unrealistically large range of electrode potentials. The dash-dotted line on the inset shows, for comparision, the simulation capacitance of a model IL where cations have the same radii as anions (0.5 nm). These data were taken from our previous work.19
24 molecular dynamic productive runs of 25 ns at constant box volume, number of particles, and temperature, preceded by 10 ns equilibration runs. In all runs the simulation temperature was kept at 450 K9 using the Berendsen thermostat.10 Each run was performed for a given charge density of the electrodes in the interval between -80 to +80 µC/cm2. The data on the ion positions were collected with 5 ps time interval. We calculated the in-plane averaged ion density profiles and potential distribution along the z-axis following a standard procedure,11-13 using a uniform grid with 0.025 nm distance between the nodes. As we have obtained about 10 nm wide electroneutrality region in the middle of the IL slab,14 the double layers at and capacitance of each electrode could be considered separately. Values of the differential capacitance were calculated by numerical differentiation of the obtained surface charge density-electrode potential dependence, with all potentials calculated relative to the electroneutral region in the middle of the slab. For more detailes on the simulation procedure, geometry of the unit simulation shell, and capacitance calculations, see ref 14. 3. Results The simulation results for the differential capacitance are shown in Figure 1 together with the theoretical predictions of EMFT3,15 for such a system.16 As one can see, these simulated curve and the EMFT curve are remarkably close to each other. The simulated capacitance curve is not Gouy-Chapman-like.
It qualitatively resembles recent experimental data on differential capacitance of the double layer in imidazolium-based ionic liquids.17 We show them for 1-methyl-3-hexylaimidazolium chloride (hmimCl) at 372 K17 in the right inset of Figure 1. Basic features of the simulated capacitance curve are 1. The capacitance has a main maximum of about 15 µ F/cm2 shifted by ∼0.3 V to the positive side of PZC as the size of cation is larger than the size of anion. As one can see, the experimental and EMFT curves have also a distinct positive shift of the main maximum. The maximum height is only slightly lower than in the drawn experimental curve and lies in the range of experimental values measured for different ILs of ∼5-60 µF/cm2.17,18 The half-width of the maximum is about 2 V, close to the EMFT results and experimental data.17 2. As in EMFT and in the experiments,17 the capacitance curve has distinct decreasing wings at large positive and negative potentials; here, the coincidence with EMFT is goodsanticipated to be so, as the asymptotic behavior of capacitance is determined by charge conservation laws (cf. ref 19). 3. The asymmetry of the curve about the PZC is also in line with the theoretical curve and experiments:17 the positive wing lies higher than the negative one. 4. Similarly to the experimental curve, the capacitance has a secondary maximum in the area of negative potentials (marked by arrow on the plot), but it is much less pronounced than in the
11870 J. Phys. Chem. B, Vol. 112, No. 38, 2008
Letters
Figure 2. Structure of electrical double layer in a model asymmetric ionic liquid near the anode (solid/blue lines) and cathode (dashed/red lines). (left column) Distributions of the volume charge density of the counter charge, displayed for the indicated absolute values of the surface charge densitites σe (µC/cm2) of the electrodes (the values shown on the right side of the figure). (right column) “Slicing” the countercharge: charge densities per unit surface area calculated for 0.1 nm thick slices. In order to assess the degree of the overscreening effect, all charge densities in this column are scaled to the corresponding absolute values of the surface charge density of the electrode. The numbers at arrows show the integrals over the first peaks and valleys at each electrode.
experimental curve. The range of potentials where the maximum emerges in the simulations is close to the experimental ones. 5. It has several secondary maxima in the area of positive potentials. They are of “academic interest” as in experiments this range is complicated by the onset of faradaic processes or electrosorption of anions. As expected, the short-range ion correlations determine the system response at potentials close to the capacitance maximum, and the EMFT theory does not exactly reproduce the simulation data, but generally, it is remarkably close to them. 4. Discussion What stands behind the dependence of differential capacitance on electrode potential? In diluted electrolytes, where the GouyChapman theory works very well, the distance of the center-ofmass of the counter charge from the electrode gets exponentially smaller with increase of electrode potential, and the capacitance rises. However, in a very concentrated ionic system, where only a limited increase of the concentration of ions in the double layer is possible, the position of the center of mass of the counter charge grows with the increase of potential (lattice saturation effect) and the capacitance falls down. For a dense, but compressible, IL the capacitance should go up with the electrode polarization as long as it draws more counterions to the electrode surface, filling the voids in the IL and thereby achieving denser packing. Once the packing is complete (all the voids filled), the capacitance will start decreasing, with counterions lining up near the electrode. There will be a maximum in between. However, in dense Coulombic systems, where the effects of ionic correlations and overscreening are important,20,21 the response cannot be rationalized solely in terms of “filling voids”
and “coming closer” of one sort of ion to the electrode. The system response is determined by the formation of standing charge density waVes, made possible by collective rearrangement of counterions and their coions. There is nothing mysterious in such patterns. While in nonpolarized system we have approximately the same amount of positive and negative ions in each lateral plane, electrode polarization stabilizes alternation of cation- or anion-rich layers, instead of alternating anions and cations in each plane. Figure 2 (left) shows the charge distribution in the double layer for different charge densities of the electrode. We see strong charge density waves, particularly for small and moderate potentials (in the range of which the main maximum of capacitance lies). Obvious trends are the following: at negative electrode potentials, the first maximum of positive charge lies father from the electrode surface than the first maximum of negative charge at positive potential, because the cations are two times larger than anions. Electrostriction is also seen due to the compressibility of the IL, which leads to the increased densities of cations and, especially, anions. To understand deeper the anatomy of the double layer, we have sliced the slab between the electrodes into the 0.1 nm layers for which we show the average charge density per unit surface area, scaled to the absolute value of the surface charge density of the electrode. The results are shown in Figure 2 (right). For small electrode charges we see that the corresponding first groups of these countercharges nearest to the electrode deliver the values that are larger than the absolute value of the charge density of the electrode. This is (over)compensated by the next layer, still excessiVely oppositely charged, and off it goes, with electroneutrality reached after a set of decaying oscillations of
Letters
J. Phys. Chem. B, Vol. 112, No. 38, 2008 11871
the counter charge. This effect is called oVescreening, for more details see refs 19 and 22. This trend breaks down at large absolute values of electrode potential where another effect gets into the play: lattice saturation. This effect is responsible for the decaying wings of the capacitance,3,23 where it approaches the asymptotic law.
C ∝ 1/√γ|U|,
U . kBT/e
(1)
with γ being the ratio of the average ionic concentration to the maximum possible local concentration. If cations and anions were of the same size, γ in MFT would not depend on electrode potential. But if the cations are larger than anions, as is usually the case for many ILs, γ must be potential-dependent. At positive electrode potentials the double layer is rich in anions, and γ is effectively smaller than in the region of compact packing of cations at negative potentials. This is the reason why the east wing lies higher than the west wing. As expected, the capacitance maximum is higher than in our previous simulations,19 where the size of anions was two times larger and equal to the radius of cations; it is also positively shiftedssee the inset on the left. The same inset shows how well the asymptotic law (1) works and demonstrates the difference in the height of the wings.24 Addition of the Stern layer contribution to capacitance, as described in ref 3 and implemented in ref 14, is instrumental for obtaining reasonable values of capacitance in the maximum. With no Stern layer “cutoff”, the maximum of the MFT capacitance is too high. But the MFT should essentially break down here anyway! Adding in series the compact layer contribution adds a “microscopic flavor” to the theory which the pure MFT inevitably misses. Such extended theory (EMFT) seems to work superbly. As in ref 19, this theory contains two adjustable parameters (see ref 14). But their values, that provide the shown good correspondence with simulation data lie exactly in the expected range. Comparing with the experiments the results of the double layer theory or simulations one should take into account that the whole capacitance curve may be additionally tilted due to the effect of asymmetry of relaxation of the surface electronic profile, widely discussed in electrochemistry of ordinary electrolytes25 but not accounted for in the classical MD simulations. That effect will provide an additional impetus for the positive derivative of the capacitance near the PZC. With the cations larger than anions, the two trends do not compete but support each other. Generally, the character of the capacitance of real ionic liquids will be affected by rearrangements of nonspherically shaped ions and more complex structure of force fields acting between them. Each ionic liquid is then expected to display “individual” capacitance features. We believe that in the future we will be able to correlate them with particular properties of ions, when comparing the data obtained for the same electrode surfaces. But this will be possible as long as we have a “point of reference”: complete understanding of the capacitance curve for model, simple ILs, like the one considered in this letter. Acknowledgment. We are thankful to Martin Bazant (MIT) for useful discussions and to Rossen Sedev (University of South Australia) for sending us a digital file with experimental data. Supporting Information Available: Description of the simulation methodology and extended mean-field theory. This material is available free of charge via the Internet at http:// pubs.acs.org.
References and Notes (1) J. Phys. Chem. B 2007, 111. (2) (a) Buzzeo, M.; Evans, R.; Compton, R. ChemPhysChem 2004, 5, 1106–1120. (b) Silvester, D.; Compton, R. Z. Phys. Chem. 2006, 220, 1247– 1274. (3) Kornyshev, A. J. Phys. Chem. B 2007, 111, 5545–5557. (4) (a) Alam, M.; Islam, M.; Okajima, T.; Ohsaka, T. J. Phys. Chem. C 2007, 111, 18326–18333. (b) Alam, M.; Islam, M.; Okajima, T.; Ohsaka, T. Electrochem. Commun. 2007, 9, 2370–2374. (5) Simultaneously with ref 3 the same expression for diffuse double layer capacitance was obtained in ref 26 for concentrated electrolytic solutions. In the context of ionic liquids, its particular case was later rederived in ref 27. In all these works, the idea of the approach was similar to those of the old, classical theory of the bulk properties of concentrated electrolytes of Bikerman,28 Dutta and Bagchi,29 and Eigen and Wicke.30 Moreover, Bazant and Kilic (private communication) have found that all the recent articles as well as their predecessors31,32 that used this kind of mean-field approach, missed a paper published in German more than 50 years ago: Freise33 derived the same formulae for the potential distribution and double layer capacitancesagain in the context of concentrated electrolytes. His work was forgotten because the predicted deviations from the Gouy-Chapman capacitance had never been observed in electrolytic solutions, as they emerge only at potentials beyond the range of ideal polarizability of the electrodes. The situation is different for ILs, where the strikingly non-Gouy-Chapman behavior should be seen well before the onset of Faraday processes. (6) Welton, T. Chem. ReV. 1999, 99, 2071–2083. (7) (a) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43–56. (b) Lindahl, E.; Hess, B.; van der Spoel, D. J. Molec. Model. 2001, 7, 306–317. (8) (a) Yeh, I.; Berkowitz, M. J. Chem. Phys. 1999, 111, 3155. (b) Yeh, I.; Berkowitz, M. J. Chem. Phys. 2000, 112, 10491. (9) We have chosen this, relatively high temperature, to facilitate the simulation convergency which is very slow for the room temperature of 300 K. We also would like to note that this temperature is close to the range of temperatures used in the recent experimental study17 of the double lazer capacitance in ILs (352-412 K). (10) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (11) Pinilla, C.; Del Po’polo, M.; Kohanoff, J.; Lynden-Bell, R. J. Phys. Chem. B 2007, 111, 4877–4884. (12) Del Po’polo, M.; Kohanoff, J.; Lynden-Bell, R.; Pinilla, C. Acc. Chem. Res. 2007, 40, 1156–64. (13) The authors of ref 11 published results of a computer simulation of an ionic liquid between charged walls, using realistic ion pair interactions potentials for dymethylimmidazolium chloride. They performed a detailed investigation of electrostatic potential relaxation in response to jumpwise charging for one value of charge density of the electrodes (σ ) 2 µC/cm2). Oscillating charge distribution profiles in a 4 nm IL-filled gap were obtained with some significant overlap due to the small width of the gap. (14) See the Supporting Information. (15) It is based on eq 20 of ref 3 with parameter γ effectively depending on electrode potential, according to eq 25 of ref 3. Next, in the extended mean field theory (EMFT), the diffuse layer contribution is combined with the Stern-like contribution. This is described in ref 3 and exploited in ref 14, but now with the width of the Stern layer changing from negative to positive potentials according to the same interpolation formula as eq 25 of ref 3 for γ. The equation for γ has no fitting parameters, because the values of γ+ and γ- are determined by the dense packing of cations and anions, respectively; the same refers to the Stern layer width, which is approximated, correspondingly, by the radii of cations and anions in the two limits. Avoiding introduction of any fudge factors, kBT was taken to be the scale of the variation of these parameters with potential. (16) The parameters used to plot the EMFT curves are the following: ˜ ) 5, ) 7, γ+ ) 0.5, γ- ) 0.07 for the diffuse layer part of the s s capacitance. The parameters of the Stern layer contribution are: d+ and dequal to 0.5 and 0.25 nm correspondingly (for a definition of these parameters and the basic formulae used see the Supporting Information). (17) Lockett, V.; Sedev, R.; Ralston, J.; Horne, M.; Rodopoulos, T. J. Phys. Chem. C 2008, 112, 7486–7495. (18) Aliaga, C.; Baldelli, S. J. Phys. Chem. B 2006, 110, 18481–18491. (19) Fedorov, M. V.; Kornyshev, A. A. Electrochim. Acta 2008, 53, 6835-6840. (20) Rovere, M.; Tosi, M. Rep. Prog. Phys. 1986, 49, 1001–1081. (21) (a) Fisher, M.; Levin, Y. Phys. ReV. Lett. 1993, 71, 3826–3829. (b) Fisher, M. J. Stat. Phys. 1994, 75, 1–36. (22) Fedorov, M. V.; Kornyshev, A. A. Mol. Phys. 2007, 105, 1–16. (23) Kornyshev, A. A.; Vorotyntsev, M. A. Electrochim. Acta 1981, 26, 303–323.
11872 J. Phys. Chem. B, Vol. 112, No. 38, 2008 (24) The height ratio between the right and left wing (or between the left wing in the asymmetric system and the left wing of the symmetric one) at the same absolute values of potential should scale, roughly, as (γ+/γ-)1/2 ) (0.5/0.07)1/2 ≈ 3. This is exactly what this inset demonstrates. (25) (a) Kornyshev, A. A. Electrochim. Acta 1989, 34, 1829–1847. (b) Kornyshev, A. A.; Spohr, E.; Vorotyntsev, M. A. Electrochemical Interfaces: At the Border line. Encyclopedia of Electrochemistry; Wiley-VCH: Weinheim, 2002; Vol. 1, pp 33-132. (26) Kilic, M.; Bazant, M.; Ajdari, A. Phys. ReV. E 2007, 75, 21502. (27) Oldham, K. B. J. Electroanal. Chem. 2008, 613, 131–138.
Letters (28) Bikerman, J. Trans. Faraday Soc. 1940, 35, 154–160. (29) Dutta, M.; Bagchi, S. Indian J. Phys. 1950, 24, 61–66. (30) Eigen, M.; Wicke, E. J. Phys. Chem. 1954, 58, 702–714. (31) Borukhov, I.; Andelman, D.; Orland, H. Phys. ReV. Lett. 1997, 79, 435–438. (32) Kralj-Iglic, V.; Iglic, A. J. Phys. II 1996, 6, 477–491. (33) Freise, V. Z. Elektrochem. 1952, 56, 822–827.
JP803440Q
