Molecular Simulations of the Electric Double Layer Structure

ARTICLE pubs.acs.org/JPCB

Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-Npropylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes Jenel Vatamanu,*,† Oleg Borodin,†,‡,§ and Grant D. Smith† †

Department of Materials Science and Engineering, University of Utah, Salt Lake City, Utah 84112, United States Wasatch Molecular Inc., 2141 St. Marys Dr., Ste 102, Salt Lake City, Utah 84108, United States § Electrochemistry Branch, Army Research Laboratory, 2800 Powder Mill Rd., Adelphi, Maryland 20783, United States ‡

bS Supporting Information ABSTRACT: Molecular dynamics simulations were performed on Nmethyl-N-propylpyrrolidinium bis(fluorosulfonyl)imide (pyr13FSI) room temperature ionic liquid (RTIL) confined between graphite electrodes as a function of applied potential at 393 and 453 K using an accurate force field developed in this work. The electric double layer (EDL) structure and differential capacitance (DC) of pyr13FSI was compared with the results of the previous study of a similar RTIL pyr13bis(trifluoromethanesulfonyl)imide (pyr13TFSI) with a significantly larger anion [Vatamanu, J.; Borodin, O.; Smith, G. D. J. Am. Chem. Soc. 2010, 132, 14825]. Intriguingly, the smaller size of the FSI anion compared to TFSI did not result in a significant increase of the DC on the positive electrode. Instead, a 30% higher DC was observed on the negative electrode for pyr13FSI compared to pyr13TFSI. The larger DC observed on the negative electrode for pyr13FSI compared to pyr13TFSI was associated with two structural features of the EDL: (a) a closer approach of FSI compared to TFSI to the electrode surface and (b) a faster rate (vs potential decrease) of anion desorption from the electrode surface for FSI compared to TFSI. Additionally, the limiting behavior of DC at large applied potentials was investigated. Finally, we show that constant potential simulations indicate time scales of hundreds of picoseconds required for electrode charge/discharge and EDL formation.

I. INTRODUCTION Understanding electrolyte interfacial properties near electrodes is important for a number of electrical energy storage systems such as Li batteries,1 supercapacitors,2 dye-synthesized solar cells,3 and fuel cells,4 as charge structuring and ion transport near interfaces influence the energy and power density of the electrochemical cells. The relationship between the electrolyte chemical structure and its ability to effectively store charge at the electrified interfaces upon applying an external potential is particularly important for the electric double layer (EDL) supercapacitors. Recently, room temperature ionic liquids (RTILs) attracted significant attention as candidate electrolytes for EDL capacitors due to low volatility,5 high thermal stability,6 wide electrochemical stability range,7 and high capacitance,8 resulting in a higher energy density. A possibility to tailor RTIL properties by choosing the appropriate cation-anion combinations makes them even more attractive for electrochemical applications, provided that the relationship between the RTIL chemical structure and its EDL properties is well understood. The shape and magnitude of the RTIL differential capacitance (DC) have recently attracted significant attention. The DC is defined as the derivative of the electrode charge per unit surface area of the electrode, Qelectrode, with respect to the applied EDL r 2011 American Chemical Society

potential, VEDL DC dQ electrode =dV EDL

ð1Þ

In 2007, Kornyshev9 pointed out that the shape of DC in RTILs is fundamentally different from dilute electrolytes and offered a theoretical treatment based on the modified Poisson-Boltzmann theory that would rationalize the observed behavior. Kornyshev9 showed how the compressibility within the double layer influences the shape of DC as a function of electrode potential. A U-shaped DC is possible when the electrolyte steric exclusion near the surface is small and therefore additional electrolyte can be packed in the surface layer with increasing potential such that the first derivative of the electrode charge with electrode potential is still increasing with potential.9,10 If the electrode surface is crowded with electrolyte, exclusion effects diminish the ability of additional ions to pack near surface and the DC decreases with increasing potential. The camelshaped DC is then, in the light of the Kornyshev model,9 a largely steric effect originating from the ion saturation (or crowding) near the surface: the increase of DC at low potentials is a signature of Received: January 5, 2011 Revised: February 4, 2011 Published: March 08, 2011 3073

dx.doi.org/10.1021/jp2001207 | J. Phys. Chem. B 2011, 115, 3073–3084

The Journal of Physical Chemistry B

ARTICLE

Table 1. Simulation Run Length and Thermodynamic and Transport Properties of ILs from MD Simulations (Density GMD, Heat of Vaporization Hvap, Ion Self-Diffusion Coefficients D, Finite Simulation Box Size Correction to Ion Self-Diffusion Coefficient ΔDFSC (See Ref 31) and Ion Conductivity λ) Using United Atom (UA) Nonpolarizable and All-Atom Polarizable (APPLE&P)31,32 Force Fieldsa T (K)

run length (ns)

FMD (Fexp) (kg m-3)

Hvap (kJ/mol)

DþMD (10-10 m2/s)

D-MD (10-10 m2/s)

ΔDFSC (10-10 m2/s)

λ (mS/cm)

expt refs

All-Atom Polarizable Force Field (APPLE&P)31,32 393

9.7

1266

1.99

2.31

0.29

42.9

333

10.1

1308

0.53

0.72

0.085

14.8 (13.4)

15

298

30.0

1339 (1339)

0.16

0.23

0.025

5.2 (6.4,8.2)

15-17

393

16.0

1239

2.28

2.87

333

4.0

1285

0.54

0.73

298

18.0

1315

0.17

0.21

143.9

United Atom (UA) Nonpolarizable Force Field

a

142.8

0.271 0.018

52.1 16 (13.4)

15

5.9 (6.4,8.2)

15-17

Experimental data are given in parentheses. The finite size correction ΔDFSC has not been included in the Dþ and D-.

loosely packed EDL that can easily accommodate additional ions as the potential increases, while the maxima and then decrease of DC at larger potential are caused by the exclusion volume effects which will eventually slow the accumulation rate versus potential of ions near surface. Recently, Kobrac11 critically analyzed Kornyshev’s treatment of DC. Lamperski showed that high densities and high temperatures favor a bell-shaped DC, while low densities and low temperatures12 favor U-shaped DC. In experimental studies, the DC was observed to be either bell-shaped,13-15 having a maxima at potentials close to potential of zero charge (PZC), or camelshaped,16,21,22 having a minimum of the U-shaped region at relatively low potentials and two maxima at higher potentials. At sufficiently large potentials the DC decreases monotonically with the applied potentials in all cases. Fedorov et al.23 and Georgi et al.24 and our previous work25 illustrated through molecular dynamics (MD) simulations that, for complex-shaped molecular ions, changes in ion orientation near the surface play an important role in explaining the shape of DC as a function of electrode potential. Specifically, weakly charged groups in the ion can play the role of “latent voids23” which, upon ion reorientation can be replaced with the highly charged groups, thus promoting U-shaped DC. Recent simulations also addressed the role of other factors, such as dispersion interactions,26 temperature dependence of DC for AgCl/water system,27 the role of electrode porosity in the storage ability of a porous electrode,28 the role of atomic polarizabilities in enhancing the crystal-like character of ordering near surface,29 and the role of electrode porosity in determining the ability of an electrode to store electrostatic energy.30 In this paper, the EDL and DC behavior for N-methyl-Npropylpyrrolidinium bis(fluorosulfonyl)imide (pyr13FSI) RTIL is studied as a function of applied potential between graphite electrodes. An extensive comparison is performed with the previously reported EDL and DC behavior of pyr13bis(trifluoromethanesulfonyl)imide (pyr13TFSI) RTIL with a significantly larger anion allowing us to focus on the influence of anion on the interfacial properties of RTILs at electrified interfaces. We also investigate the limiting behavior of DC versus the electrode potential U, and check if the U-0.5 asymptotic decay as predicted by Kornyshev9 model holds for pyr13FSI. In the employed simulation setup, we probe relaxation times of EDL of hundreds of picoseconds. The paper is organized as follows: in section II we briefly describe development and validation of an efficient united atom

force field for pyr13FSI that will be used in the simulations of pyr13FSI sandwiched between graphite as described in section III together with the simulation methodology. In section IV, behavior of DC is presented as a function of applied voltage, while the EDL structure and its correlations with DC are analyzed in sections V and VI and are related to the DC behavior. The asymptotic behavior of the DC at high voltages is discussed in section VII. The kinetics of electrode discharging is presented in section VIII. Finally, concluding remarks will be given in section IX.

II. IONIC LIQUID FORCE FIELD DEVELOPMENT In order to reduce computation costs, an accurate nonpolarizable united atom (UA) force field for N-methyl-N-propylpyrrolidinium (pyr13þ) bis(fluorosulfonyl)imide (FSI-) was developed and validated as described in detail in the Supporting Information. Here, we present a brief summary of the ability of MD simulations using this UA force field to describe heat of vaporization, density, and ion transport in pyr13FSI. Table 1 indicates good agreement between MD simulation results using the UA force field with experimental data and previous predictions made using the more computationally expensive APPLE&P many-body polarizable force field.31,32 This agreement is superior to the previous predictions using a nonpolarizable all-atom version of the APPLE&P force field.33 Specifically, a nonpolarizable all-atom force field for pyr13FSI created by turning off atomic polarizability predicted a significantly higher heat of vaporization of 177.8 kJ/mol compared to the value of 143.9 kJ/mol value obtained for APPLE&P.33 In contrast to these results the current UA force field predicts a heat of vaporization only 1 kJ/mol different from APPLE&P. Ion diffusion was also a factor of 2.6-3.0 slower at 393 K in the nonpolarizable all atom version of the pyr13FSI force field, while the current UA force field predicts ion transport within 25% of APPLE&P values. We conclude that the UA nonpolarizable force field for pyr13FSI developed in this work provides an adequate description of pyr13FSI bulk properties. III. SIMULATION METHODOLOGY In this paper, we present a molecular dynamics (MD) simulation study of the EDL structure and differential capacitance for pyr13FSI sandwiched between graphite electrodes. The chemical structure of the ions indicating the atom labels utilized throughout this paper and an image snapshot of our simulation setup is shown in Figure 1. The simulation cell consisted of 190 3074

dx.doi.org/10.1021/jp2001207 |J. Phys. Chem. B 2011, 115, 3073–3084

The Journal of Physical Chemistry B ion pairs and two layers of graphite oriented with the basal face toward electrolyte. Each graphite electrode consisted of 480 atoms. The cross-sectional area of the system was 25.614 24.647 Å2. The asymmetry direction, which is defined as the direction perpendicular to the graphite surface, is denoted as the z-axis. The size of the simulation cell was chosen to be large enough in order to have no ion density and potential oscillations at the center of the cell as shown in Figure 2. The distance between the electrodes for pyr13FSI was 125.8 Å, which ensures sufficient bulk phase between the two electrodes. In order to perform head to head comparison with our previous work

Figure 1. A schematic representation of (a) TFSI anion, (b) FSI anion, and (c) pyr13 cation. The labels are the atom notation used in the paper. The arrows indicate the molecular axes used to compute the orientation of the ions relative to electrode surface. The bottom figure represents a snapshot of a typical simulation setup for a 2D-geometry system. The electrode atoms are represented as gray spheres.

ARTICLE

pyr13TFSI-graphite,25 the pyr13FSI-graphite simulations were also performed at 393 and 453 K. Dimensions of the simulation cell were adjusted to yield density in the middle of the simulation cell equal to the bulk density of pyr13FSI at 1 atm. The bulk (3D) simulations for force field development and validation were run at 298, 333, and 393 K, while 2D simulations of the confined IL were performed at 393 and 453 K. The Poisson potential, φ(z), was obtained by integrating the Poisson equation d2 φðzÞ=dz2 ¼ -FðzÞ=ε0 for the measured charge profile along the z-direction, F(z). The Poisson potential was found to be flat in the middle of the simulation cell, indicating zero electric field in the middle of the cell, and hence a complete screening of the applied field on electrode by EDL. In order to compute the DC at any given temperature, we generated 32-40 ns trajectories at potential differences between electrodes, ΔU, ranging from to 7 V. The simulations were performed under constant potential difference between the two electrodes, in a similar setup as described in ref 34 and in our previous paper ref 35. The electrode charges were computed iteratively, minimizing the total electrostatic energy subject to the constraint that the electrostatic potential on the positions of the electrode atoms is the same and equal with the dialed potential. The electrode charges were represented as Gaussian distributed36 having widths of 0.5 Å.37 The simulation setup consisted in two graphite electrodes with the basal face oriented toward the electrolyte, and a RTIL electrolyte inserted between electrodes. Note that the addition of the second layer of graphite is a small refinement; a single graphite layer should be sufficient to reproduce the results presented here. The position of the electrode atoms was frozen during simulations. The short-range nonbonded forces (i.e., van der Waals and real part of Ewald) were evaluated within a spherical cutoff of 10

Figure 2. Ion center of mass density profiles at 393 K for pyr13FSI and pyr13TFSI at 393 K (a) at PZC and (b) at ΔU = 5.2 V applied potential between electrodes. 3075


The Journal of Physical Chemistry B Å, while the long-range electrostatics were computed with smooth particle mesh Ewald adapted for 2D geometry.38 The equation of motions were integrated with a symplectic multipletime-step algorithm39 over time resolutions of 0.166 fs for bonds, angles, and out-of-plane deformations, 0.833 fs for dihedrals and nonbonded short-range interactions within 7.5 Å cutoff, and 5.0 fs for complete update of long-range forces. The electrode charges and the measured quantities were updated every 0.1 ps. The temperature was controlled with a Nose-Hoover thermostat40 having a coupling constant of 0.5 ps.

IV. DIFFERENTIAL CAPACITANCE AND POTENTIAL OF ZERO CHARGE We begin our analysis with discussing the PZC and DC of EDLs on the positive and negative electrodes. As PZC is nonzero, we set our electrode potential scale (URPZC) as the drop of

Figure 3. Differential capacitance versus the electrode potential relative to PZC for pyr13FSI and pyr13TFSI.25

ARTICLE

potential within the EDL relative to PZC: URPZC ¼ UEDL - PZC ¼ φelectrode - φbulk - PZC Note that throughout this paper we will use a simplified notation URPZC designating the number of volts of a potential scale from which PZC was subtracted out. The PZC for the pyr13FSI system was found to be -0.560, -0.574, and -0.554 V for systems at 363, 393, and 453 K, respectively. The PZC was determined as the difference between the Poisson potential on electrode and in bulk for the simulations at 0 V potential differences between the two electrodes. The negative PZC indicates a larger affinity of the surface toward anions, due to the presence of exposed electronegative oxygen and fluorine atoms in the anion next to the electrode as discussed below. For pyr13TFSI, we found the PZC to be -0.34 and -0.41 V for systems at temperatures of 393 and 453 K.25 The larger PZC for pyr13FSI versus pyr13TFSI is due to the higher charge of the anion atoms situated immediately next to graphite surface: in the graphite|pyr13FSI interface, oxygen atoms are located next to graphite, while in the graphite|pyr13TFSI interface much weaker charged fluorine atoms are located the closest to the graphite surface. The DC was evaluated via eq 1 using a local fit of the chargepotential dependence shown in the Supporting Information as follows: for each set of 11 (or 13) consecutive points EDL potential-electrode charge, a parabola was fitted. Then the DC was computed in the midpoint (the sixth point) as the numerical derivative of the locally fitted parabola. The procedure was repeated from all 63 pairs of data electrode charge-EDL potential. Shown in Figure 3 is DC versus EDL potential relative to the potential of zero charge (VRPZC), for temperatures of 393 and 453 K. An overall camel-shaped trend of DC is observed at 393 K, with maxima in DC at around -1.5 and þ1.7 VRPZC. Compared with pyr13TFSI, the DC for pyr13FSI is larger at the

Figure 4. Density profiles of various atoms. These densities are normalized by the number of atoms of the same type in the molecule. For example, plotted for O groups is the total density divided by 4 atoms of O in FSI molecule. The potentials are relative to PZC. Represented in the X axis are the distances from the electrode in Å, in the Y axis are the EDL potential relative to PZC, and in the Z axis are the actual values of the densities. 3076


The Journal of Physical Chemistry B negative electrode and has the position of its minimum shifted toward positive EDL potentials. The U-shaped region almost disappears at 453 K such that in the range -3.5 to 0 VRPZC the DC is bell-shaped (with its maximum at -1.5 V), followed by an almost constant plateau in the potential range 0 to þ2.5 VRPZC. Beyond the maxima in DC, we observe a steep decrease, of about 2-3 μF/cm2 per volt, at both positive and negative electrodes. In the next two sections we describe EDL structural changes as a function of applied voltage and relate them to the observed DC behavior.

V. ELECTRIC DOUBLE LAYER STRUCTURE A comprehensive characterization of the electrode-electrolyte is provided by the density profiles along the asymmetry direction (z-axis) and the distributions of the ionic orientation relative to surface of the innermost ions. V.1. Ion Density Profiles. The profiles of density, charge density, and Poisson potential in the z-direction (perpendicular to surface) show multilayer ordering near the electrode extending up to 20-30 Å at PZC and up to 40 Å for the cell with 5.2 V applied between electrodes. The density profiles along the entire simulation box are shown in Figure 2 while the EDL region is shown in Figure 4. Behavior of N(pyr13) and N(FSI) in Figure 4 is similar to the behavior of the center of center of mass density profiles, and therefore the center of mass density profiles are not shown in Figure 4. The center of mass density profiles near the electrodes exhibit an increase of counterions and a sharpening of their distribution, as the electrode potential increases while co-ions move out of the innermost layer of EDL. The first peak of the pyr13 density profile moves toward the electrode by about 0.5 Å as the potential is increased from PZC to large negative potential (≈-3.5 V for EDL on the negative electrode). The position of the first peak of the FSI density does not significantly change with the positive potential increase from PZC to large positive potential (≈þ3.5 VRPZC for EDL on the positive electrode) while its magnitude significantly increases. On the basis of these density distributions, we define the inner and intermediate layers of EDL. The inner layer extends from 0 to 5.3 Å from the graphite surface and encompasses the center of mass of counterions adsorbed on the graphite at high potentials as shown in Figure 4. The intermediate layer encompasses the electrolyte distribution in the range 5.3-7.5 Å from surface and captures most of the second maximum in density profile distributions. The co-ions move from the inner layer into the intermediate layer and beyond with increasing EDL potential. At high negative and positive voltages the center of mass of counterions are located primarily in the inner layer, while the center of mass of co-ions are located in the intermediate layer and beyond. Figure 4 shows additional details for the atomic density profiles in the EDL region as a function of distance from the electrode surface for the entire range of electrode potentials investigated. The two maxima in the density profiles of CN (pyr13) located at 4 Å and 6 Å from surface observed at potentials more negative than -1 VRPZC are consistent with perpendicular to the surface orientations of the 5-ring cycle (see next). The O(FSI) and F(FSI) density profiles have two well-defined peaks near graphite electrode at ∼3.3 Å and 5.3 Å that become very sharp at high positive potentials. This behavior is explained by the parallel to surface orientation to surface of FSI such that one fluorine atom points toward surface and the other away from surface, or arrangements where both fluorine atoms point away

ARTICLE

from the surface and four oxygen atoms are bonded to surface (see image snapshots in Figure S2 in the Supporting Information). The density profiles of Figure 4a indicates that the first peak of the CH3 group of pyr13 shows an increase of the probability of finding CH3 groups of pyr13 near the surface as the potential becomes more negative. Interestingly, even at large positive potentials (≈þ3.5 VRPZC) the terminal CH3 groups can penetrate near the surface, although the N(pyr13) atom shown in Figure 4b and therefore the mass center of the cation are repelled from the positive surface. A more detailed analysis of the anion orientation is presented in the next section. While the largest variation in structure occurs within 15 Å from the surface, three to four additional layers in density oscillations were observed extending over 35-40 Å from the electrode (Figure 2). At high potential, the minima in cation density profiles coincide with the maxima at the anion density profiles, indicative of ion multilayer formation. In spite of the complex molecular structure of pyr13FSI, the oscillations in EDL structure beyond 15 Å from surface can be described by a fairly simple function and are similar to those for simple (atomic) ionic liquids:35 Δz Δz þθ ð2Þ FðΔzÞ ¼ K exp sin 2π λ β where Δz is the distance from the electrode, λ is the screening length which indicate penetrability of the electric field in electrolyte, β shows the periodicity of the charge oscillations and it is related to the size of the ions(see also the inset of Figure 2), θ is a phase factor, and K is a proportionality constant dependent on the magnitude of the applied potential on electrode. Within the errors of our fit, we find that the screening length λ is independent of the applied field, and it decreases with temperature increase. For both systems the screening length λ is 12.0 Å at 393 K and 11.1 Å at 453 K. A decrease of screening length with temperature is expected because ordering is favored by lower temperature. The periodicity of oscillations β is very similar for FSI and TFSI systems, specifically of 7.0 Å and 7.3 Å, respectively, at 393 K. Slightly larger β for TFSI system is expected due to larger molecular volume of TFSI. V.2. Distributions of Ionic Orientations Relative to Surface. Further details of the EDL structure are provided by the distributions of the ionic orientations relative to electrode surface. These distributions show a consistent picture with density profiles disused above. The distribution of the ionic orientations relative to the z-axis (perpendicular to the electrode surface) was analyzed for ions within 7.5 Å from the electrode surface. For comparison purposes with our previous work, we used the same definition of the innermost and intermediate layers as described in the previous section. For FSI, the orientation was shown for the elongated axes of the molecule defined as the vector connecting the two S atoms (see the arrow from Figure 1a). For pyr13 we used the “ring vector” connecting the middle of CU-CU bond with N atom (see the arrow in Figure 1b). These distributions are weighted by the ion population in the layer and computed by binning the orientations in equal intervals in the cosine of the angle between the molecular axis and the axis normal to surface. We verified that for bulk electrolyte, i.e., ions located sufficiently far from the surface, these probabilities show no preference for any angle of orientation, as expected. The summation of the probabilities over layers and the ionic species is 100%. 3077



ARTICLE

Figure 5. Distribution of the orientation probability with OZ direction of the elongated axes of pyr13 (right panels) and FSI (left panels) at PZC (upper panels) and large positive or negative electrode potential (lower panels). The EDL potentials (relative to PZC) are indicated in legends. The black and red lines are results for the innermost and intermediate layers. The solid symbols are results for pyr13/FSI system. For comparison we are shown results for the previous pyr13/TFSI system, (ref 25) represented as open circles and renormalized over population within 7.5 Å from electrode. The distributions of FSI/TFSI from 90 to 180° are symmetric with distributions in the range 0-90°, hence only one quadrant (0-90°) is shown in the left panels. The numbers on the OX axes are the angles between the elongated ionic axes and the direction normal to surface. The numbers on the OY axes are probabilities.

In order to illustrate the dependence of ionic orientations on the applied potential, we present in Figures 5 and 6 the distributions of electrolyte orientations relative to the z-axis, at different potentials: at PZC and at large potentials of ≈(3.5 VRPZC. The representative snapshots of ion orientations in the inner and intermediate layers of EDL are shown in Figure 7 in order to visually aid interpretation of ion orientations shown in Figures 5 and 6. There are several similarities in orientation distributions between pyr13FSI and pyr13TFSI EDLs: (i) these orientations are most structured in the innermost layer, (ii) for both systems we observe an overall increase of the parallel to surface orientations of the innermost ions, as the potentials increases, (iii) near PZC, the orientations of the ring vector of pyr13 ions located in the intermediate layer are very similar for both systems (see Figure 5b), and (iv) the intermediate layer pyr13 contains, at PZC, a large population of perpendicular to surface orientations were the alkyl tail points toward the surface; note that such orientations were also observed in a series of other experimental observations41,42 and recent simulations on a similar RTIL.43 There are, however, several surprisingly large differences in orientations of FSI and TFSI ions near electrode, mainly due to the fact that the smaller molecular volume of FSI allows less constrained packing than the more voluminous TFSI ion containing large weakly charged CF3- groups, at otherwise similar surface densities. These differences are summarized as follows:

The FSI orientations beyond 5.3 Å from surface are essentially random. This is expected because FSI is less elongated than TFSI. (ii) The distributions of the innermost pyr13 orientations on the negative surface are sharper for pyr13FSI than for pyr13TFSI (see Figures 5b,c and 6). (iii) At PZC, the distributions of orientations are overall broader for the FSI system. Specifically, at PZC, the distributions of orientations of the innermost FSI anions have two maxima: one corresponding to a tilted orientation at about 35° relative to the surface and the second parallel to surface, while the distributions of orientations for the TFSI anion only have one single maximum (see Figure 5 and ref 25). (iv) At PZC, the innermost pyr13 has, for the FSI system, two preferred orientations relative to surface: perpendicular (with the alkyl tail pointing away from the surface) and a tilted by 25° orientation, as shown in Figure 5b. In the case of the TFSI system, the innermost pyr13 was tilted by about 35° relative to the surface.25 A consequence of such difference in orientations is a significantly larger pyr13 population within 5.3 Å from surface (our defined innermost layer) for the FSIpyr13 system. (v) At large negative potentials (i.e., below -3.4 VRPZC) the co-ion is essentially completely repelled from the negative surface, and thus we expected the arrangement of the innermost pyr13 to reflect the inherent packing ability of pyr13 and be largely independent of the nature of the coion. However, as shown in Figures 5d and 6, the

(i) The distributions of the innermost FSI orientations on the positive surface are broader than for TFSI (Figure 5a,c). 3078



ARTICLE

Figure 6. A comparison of the distributions of angle between 5-ring cycle and electrode surface for the TFSI/pyr13 system and FSI/pyr13 surface for the pyr13 ions located in the innermost layer.

orientations of the 5-member ring of pyr13 in the innermost layer relative to the electrode surface are surprisingly different for these two systems at -3.45 V. For pyr13TFSI the distribution of the 5-member ring orientations relative to surface is much broader than for pyr13FSI at similar potential, having three maxima at tilts of 20°, 80°, 120° with respect to the plane surface, while in the case of pyr13FSI, the pyr13 is more tightly packed and the 5-member ring of pyr13 is oriented perpendicular to surface, as shown in Figure 6.

VI. CORRELATION BETWEEN EDL STRUCTURE AND DIFFERENTIAL CAPACITANCE MD simulations give us a unique opportunity to explore the similarity and differences of DC behavior of pyr13FSI and pyr13TFSI in terms of changes of EDL structure. We attribute the decrease of DC with increasing the absolute part of the potential beyond 2.5 VRPZC to the formation of well-defined double layers (as shown in the Figure 2), and counter-ion saturation in the innermost layers. Both RTILs also show U-shaped DCs at relatively low potentials. According to Poisson-Boltzmann theory,9 a U-shaped DC is expected at low temperatures12 and low densities.9 We estimate that our system is close to the threshold of transition from U-shaped to bellshaped DC with the cation-anion binding energy around 80 kcal/mol, which is around 100 kBT at 393 K. Presence of both weakly and strongly charge groups together with the changes of asymmetric ion orientation with potential in the innermost layer would favor the U-shaped DC, if such rotations remove weakly charged groups from surface and replace them with stronger charger ones.24,25 However, as the FSI ion is more spherical than TFSI, its changes in orientations with potential within EDL would be less important in determining the shape of DC. The asymmetry in size between FSI and pyr13 could be expected to favor U-shaped DC, if the replacement of one large ion (pyr13) from surface generate sufficient empty sites on surface to allow large rates (vs. potential) of smaller ions (FSI) accumulation; it is, however, unclear if FSI and pyr13 are sufficiently different in size to generate U-shaped DC from pure phenomenological considerations.9 Understanding the difference in DC in terms of differences in EDL structure for the two systems is less straightforward. Here we consider the following factors important in determining the shape of DC to explain the observed difference: (a) the extent

Figure 7. Representative snapshots and schematic depictions of the interfacial ion orientations for pyr13FSI RTIL and pyr13TFSI at 393 K for five EDL potentials. In these snapshots, we selected ions whose centers of mass were located within 8 Å from the electrode surface. See Figure 1 for atom labels. The electrode atoms are in gold color. The green vertical lines schematically represent the electrode surface.

(versus potential) of desorbing of co-ions and adsorbing counterions with potential increase; (b) how compact the ions can pack near the electrode surface, and (c) how close to the surface the counterions can get, specifically, the closer the screening charge is located to the electrode surface, the larger DC is expected to be. The locations of the pyr13 and anion center of mass relative to the graphite electrode are compared in Figure 8 for pyr13FSI and pyr13TFSI at 393 K. Both cation and anion approach graphite closer in pyr13FSI compared to pyr13TFSI RTILs at PZC, consistent with the greater DC for the former at PZC (see Figure 3). The pyr13 cation is located 2 Å closer to the surface, contributing to the larger DC values observed on the negative plate for pyr13FSI system than pyr13TFSI one. Next, we attempt to understand the contribution of the individual groups located within the EDL to charge screening that is connected to DC. We analyzed the derivatives with respect to potential of the average densities within various layers near electrode. The average densities in the innermost layer of O and N(pyr13) are shown in Figure 9. These average densities were calculated by dividing the average number of ions (or atoms) in a given layer to the volume of the layer. The data available for pyr13FSI allowed us to extract the numerical derivatives with respect to EDL potential of the layer density for each individual 3079



ARTICLE

Figure 8. Comparison of ion densities near electrodes for pyr13FSI and pyr13TFSI RTILs at PZC and at negative and positive potentials. All potentials shown in the legend are relative to PZC.

atom, and the contribution to total charge from each atom and layer. Before numerical differentiation, we equally spaced the data on an equidistant potential scale by linear interpolating the crude MD data, and further smooth it with 3-point adjacent averaging procedure two times. Shown in Figure 10 are the derivatives with respect to the electrode potential of the charge (a) within 5.3 Å and (b) 7.5 Å from electrode surface and the contributions to total charge from O atoms, F þ O atoms, and F þ O þ N(pyr13) atoms. The derivative of the total charge within both distances resemble the main features of DC, specifically, a minimum near PZC, larger derivatives on the negative electrode for FSI, and a drop in these derivatives beyond surface saturation at (2.3 to (2.5 VRPZC. This suggests that, in spite of complex multilayer EDL structure, DC shape is largely determined by the electrolyte closest to the surface, i.e., within less than 7.5 Å from the electrode. The derivatives of the O charge show a minimum near PZC, which reflects the low tendency of O to leave the surface with decreasing potential for the range -0.1 to þ0.3 VRPZC (see also the small slope in the O(FSI) density near PZC in Figure 9). Adding to the O-charge derivatives the contributions from other atoms (such as F and N(pyr13) shown in Figure 10) results in curves which qualitatively resemble the shape of both O-charge derivative and DC. Therefore, we can trace the U-shaped character of DC in the minimum of the innermost O density derivative (Figure 10). This suggests that a stronger adsorption to the surface of the innermost FSI via O atoms observed at low potential would be an important factor in lowering the DC near the PZC and, therefore, determining the U-shaped DC at low potentials. Higher concentrations of innermost (within 5.3 Å from surface) pyr13 and faster rates with negative potential increase of the negatively charged O(FSI) removal from surface (see Figure 9) are likely the main factors determining the larger DC for pyr13FSI on the potential range ≈-2 VRPZC to þ0.5 VRPZC. Note that larger concentrations of the innermost pyr13 for pyr13FSI are consistent with less random (i.e., sharper) distributions in orientations of the pyr13 shown in Figures 5 and 6. Larger DC near PZC for pyr13FSI system could also be

Figure 9. Average densities for N of pyr13þ and O of anions for pyr13FSI and pyr13TFSI RTILs, for atoms located within 5.3 Å from the electrode surface.

correlated with differences in orientations near PZC of the innermost anions; specifically, more parallel to the surface orientations for FSI than TFSI will place more O groups on surface for FSI anion. The more compact packing of pyr13 near the negative surface and its higher density in the innermost layer observed for the pyr13FSI system would explain why the negative electrode saturates faster for the FSI system than for the TFSI one. At the positive electrode, the FSI system requires larger positive potentials than pyr13TFSI to saturate. The smaller footprint of FSI compared to TFSI has the net effect of allowing co-ion groups such as N, CU, and CN to penetrate and fill the positive electrode surface (see Figure 9 and Supporting Information), extending the potential range at which the positive surface is crowded only with FSI counterion. We expected larger DCs on the positive plate for FSI because, due to its smaller volume, it could pack more densely. However, as shown in Figure 3, the DC on the positive plane in the potential range þ0.5 to þ1.8 VRPZC are very similar for the two systems. This can be explained as follows: the increase in negative charge in the innermost layer for TFSI system is mainly due to 3080



Figure 10. Derivatives of total charge (blue dot-dashed line), and the corresponding contributions to charge from O(FSI) (black solid line), O þ F(FSI) (red dashed line), and O þ F þ N(pyr13) (green dotted line), calculated counting the electrolyte over distances of (a) 5.3 Å and (b) 7.5 Å from the electrode surface, for the pyr13FSI system. For easier comparison with DC, we rescaled these derivatives in units of μF/cm2 and set them to have the same sign as the derivatives of the electrode charge.

ARTICLE

Figure 11. Dependence of electrode charge versus the square root electrode potential, for potentials larger than 4 V.

the increase of TFSI concentration at the surface, while for FSI it is due to a combination of slower negative charge adsorption and faster positive charge desorption with the net effect of generating comparable DCs for the two systems on the positive plate.

VII. ASYMPTOTIC BEHAVIOR OF THE DIFFERENTIAL CAPACITANCE Upon charging of the electrode at sufficiently large potential, the electrode surface becomes saturated with counterions and a growth of the effective width of the double layer is expected9 rather than additional crowding on the innermost layer. Accounting for this effect, the Kornyshev model9 predicts, in the saturation regime, a decrease of DC with potential as DC µ U-0.5 and correspondingly, an electrode charge increase as µU0.5. This behavior was explained by Kornyshev from the principle of charge conservation on a saturated double layer which can grow in width as Leff µ U0.5. We expect this behavior to be general and system independent as the assumptions made (lattice saturation and charge conservation) are fundamental. Such behavior was verified for simple models in the original Kornyshev paper, in the subsequent work of Fedorov et al.44,45 for electrolytes consisting in different sizes of ions, and in the recent work of Georgi et al.23,24 for coarse-grained molecular ILs. In contrast, we have shown in recent work25 that the DC for pyr13TFSI decreases faster than U-0.5 over a rather large potential window above the maxima in DC, i.e., in the range þ1.7 to þ3.5 VRPZC. A similar faster drop in DC beyond the potential of DC maximum was also reported by Lauw et al.46 for dendrite-shaped ILs. In order to further investigate the behavior of DC at high electrode potentials, we ran several simulations of pyr13FSI in the range of 3-50 V. For practical purposes, we emphasize the rather academic interest of the simulations at such large potentials. Also, we remind that, as the DC is the derivative electrode charge versus potential, steeper than U-0.5 decrease of DC corresponds to a slower than U0.5 increase of the electrode charge. These simulations were performed in a slightly different manner than simulations at lower potentials. Specifically, we distributed the charge equally over the surface layer of basal graphite and generated two sets of 6 ns simulations. The Poisson potential drop was computed from the charge distribution within the EDL for each electrode. To ensure that the electrodes are well separated, we increased the distance between electrodes to about

Figure 12. Profiles of atomic mass densities at large potentials. The inset indicates the preponderant orientation relative to surface of pyr13 located in the first sharp peak (i.e., innermost layer). The surface atoms are in gold color.

390 Å, and correspondingly, the amount of bulk fluid between electrodes. Shown in Figure 11 is the absolute value of the electrode charge versus square root of the absolute value of the EDL potentials in the range of 4-49 V. The electrode charge clearly exhibits a linear dependence on the square root of potential in the potential range between 15 and 49 V, thereby following the expected DC µ U-0.5 behavior for the highly charged surfaces. However, in agreement with prior results25 for pyr13TFSI over a smaller potential range, the electrode charge for the 4-10 V potential window exhibits a slower than U0.5 increase, and therefore the DC decreases faster than U-0.5 on this potential range. Density profiles and image snapshots of the electrode/electrolyte interface provide an explanation for the faster than U-0.5 decrease of DC. In Figure 12, we show the profiles of 3081


The Journal of Physical Chemistry B center-of-mass molecular density near the negative electrode. While at and below -10.82 V there are still alternating cation and anion layers near the electrode, at negative potentials above -19.45 V the width of the innermost layer counterions grow (e.g., compare Figure 12, b and c). In other words, we observe overscreening up to surprisingly large potentials (≈ (15 V), and overcrowding (i.e., a growth of the innermost layer thickness) above (15 V. It is interesting to note the surprisingly large potential (about 15 V) required for transition from overscreening to overcrowding observed for our system. For comparison, Bazant et al.47 found for simple spherical models much smaller potentials, (2.6 V) for this type of transition. The large binding anion-cation energies, of about 100kT in the gas phase, would favor overscreening for the molecularly complex RTIL. The transition from overscreening to overcrowding is associated with other structural changes in EDL described below. For conformationally flexible ions such as pyr13, conformational changes with increasing electrode potential can be expected, and some conformers preferred near the electrode at large voltages could require a change in packing in innermost layers, resulting in a change in EDL structure rather than the asymptotic growth of the EDL width. This effect could generate a slower than asymptotic growth of the electrode charge, if such changes within EDL require high activation energies, and occur slowly with potential increase. For example, in the potential range beyond maxima in DC and up to -10.82 V, the innermost pyr13 cations prefer to align with their 5-member rings perpendicular to surface plane, and the alkyl tail roughly parallel to surface plane (Figure 12a). Below about -15 V, the 5-member ring prefers to reorient such that the ring becomes parallel to the plane, while the alkyl tail is repelled from the surface (see Figure 12b). Note than Georgi et al.24 also predicted neutral (or weakly charged groups) pointing away from the surface at large potentials, but in their case a rotation of the ion was sufficient to reorder the tail, and much lower potentials were required (about 3 V). In our case, the need to assume a stretched tail conformer (as shown in Figure 12c) may explain the much larger potentials at which neutral tails are repealed from surface. At potentials more negative than -20 V, the innermost pyr13 is predominantly in a highly stretched configuration and this particular configuration is maintained in the innermost layer over the entire potential range of 20-49 V which exhibits the asymptotic behavior DC µ U-0.5. To conclude, the faster than U-0.5 intermediate decay in DC observed here corresponds to the potential window between maxima in DC and the transition between overscreening to overcrowding. At sufficiently large potentials, when the structure/conformation of the adsorbed ions remains unchanged at potential increase and the inner layer becomes overcrowded and growing in width, the DC µ U-0.5 regime is reached for the complex RTIL studied here, a fact strongly supporting the underlining fundamental assumptions behind the simple EDL models that predict this asymptotic behavior.9 Finally, we note that, in a recent paper published while ours was in revision, Bazant et al.47 identified, with a modified Landau-Ginzburg-type continuum theory, an intermediate asymptotic decay of DC as U-1/4. As for the particular systems discussed here, we only identify steeper than U-1/2 decay of DC; at this point the issue of intermediate decay of DC would still seem perhaps controversial. Further simulations on various RTIL are clearly warranted to clarify this aspect. For example, RTIL with longer alkyl tail would appear to be a good candidate to test Bazant’s slow intermediate decay.

ARTICLE

VIII. KINETICS OF ELECTRODE CHARGE AND ELECTROLYTE-INDUCED FIELD RELAXATION With a suitable methodology, MD simulations can provide valuable atomistic insights into the mechanism of EDL formation. In this section we will show that our constant potential simulations reproduce relaxation time scales of the EDL formation in good agreement with the experiment. We show that the discharging of a highly conductive electrode is limited by the diffusion processes near the electrode and not by fast (i.e., optical) vibration modes, and we probe a time scale required for electrode charging or discharging on the order of hundreds of picoseconds. The dynamic studies reported in this section were performed with two different methods: (i) constant charge simulations, where the electrode charges were suddenly dropped to zero and then maintained constant, and (ii) constant potential simulations, where the potential between the electrodes is suddenly dropped to zero and then maintained at ΔU = 0 V. In the second case, the electrode charges fluctuate such that the electrostatic potential on each electrode atom position is zero and therefore the charge distribution on the electrode would be expected to mimic the electrostatic distribution on a highly conductive electrode. As we will show, the constant charge setup predicts a relaxation of large induced electrostatic field over time scales of picoseconds while the constant potential setup predicts a time scale of the electrode charging/discharging and EDL formation of hundreds of picoseconds. In the constant electrode charge simulations, we discharged an initially charged system of two electrodes by suddenly setting the electrode charges to zero, and allowed the system to evolve without any constraints on the electrode potential. In the initially charged system, the electrolyte was subject to a large external field given by the potential drop between the two electrodes in vacuum. The electrolyte placed between the electrodes responds to such external field generating an induced field in the electrolyte which cancels out the external field such that the total field in bulk, for an equilibrated system, is zero along our z-direction. When the charges on the electrode are suddenly removed, the external field is also removed and we will observe only the induced field in the electrolyte and its evolution in time. Therefore, in the constant charge setup we study the relaxation of a large electrostatic field suddenly induced in the probe. A similar setup was previously employed by Pinilla et al.48 for dimethylimidazolium chloride, and Lanning et al.49 for molten KCl. The dependence of the induced field versus time after external field was removed was computed by integrating the electrolyte charge distribution for instantaneous configurations and averaging over 10 independent runs. In agreement with the results of others,48-50 we found that the relaxation of these large electrostatic fields can occur very fast, over several picoseconds, in two steps: a fast initial decay over 1 ps, followed by a slower second step relaxations occurring over about 3 ps. Clearly, over such small time scales, the ion diffusion has little impact. Thus, the large fields suddenly induced in the electrolyte will relax via optical mode vibrations, as proposed by Pinilla48 and Lanning.49 Consistent with the results of Pinalla,48 current densities along the z-direction defined as jðtÞ ¼ Σþ qþ ðdzþ =dtÞ þ Σ- q- ðdz- =dtÞ where q is the ion charge and dz/dt is the z-component of ion velocity, show for our system the same sign for both ionic species, 3082



ARTICLE

Figure 13. (a) Relaxation of the electric field induced in the electrolyte. The inset represents the variation in time of the current densities of each individual ion and of the total current. (b) Relaxation of the electrode charge after the potential difference between electrodes is suddenly changed and then maintained at zero.

indicating field relaxation via fast vibration modes that involve motion in opposite directions of the ionic species of opposite sign charge (see the inset of Figure 13a). While the constant charge setup verifies that large induced electric fields can relax within picoseconds,48 a different simulation approach should be used in order to capture the realistic time scale of EDL formation. Indeed, experiments indicate time scales of hundreds of picoseconds required for EDL formation51 as typical for RTILs. We will show that the constant potential simulation setup, where the electrode charges are updated via image charge method as described in ref 34, predicts realistic time scale for electrode charge relaxation and EDL formation. We started from systems charged at ΔU = 3, 4.2, and 6 V potential difference between the electrodes, and suddenly dropped the potential difference between the electrode at ΔU = 0 V. During discharge, we imposed constant and controlled potential difference between electrodes (of ΔU = 0 V) and allowed the electrode charges to fluctuate according to the imposed potential and the electrostatic environment. In these simulations, we used a time step of 3 fs and updated the electrode charges every time step. Shown in Figure 13b is the variation of the electrode charge in time, after the potential was suddenly dropped and then maintained at 0 V. Clearly, the time scale over which the electrode discharges, i.e., the electrode charges change to values corresponding to the new imposed potential, is of several hundreds of picoseconds. As shown in Figure 13b, there is an initial fast decay in electrode charges occurring in 200 ps, followed by a slower decay over about 1000 ps. The total relaxation time of the electrode charges, of hundreds of picoseconds, is long enough to allow the rotational and translational diffusion to reorder the ions within the EDL according to the new imposed potential. We verified from image snapshots that the innermost ions were indeed reoriented after about 1 ns, according to the expected structure at PZC. Under the constant potential setup we find (as expected) for highly conductive electrodes that the electrode discharge induced upon a sudden change in the electrode potential is limited by the diffusion processes that reorganize EDL.

IX. CONCLUSIONS In this paper, we present MD simulations of the EDL structure and DC for pyr13FSI RTIL near the basal face of graphite, as a function of temperature and the applied potential. A camel-like shape of DC was observed for pyr13FSI at 393 K. Larger values of DC observed on the negative electrode for pyr13FSI compared to pyr13TFSI are in agreement with recent experiments reporting larger DC on the negative electrode for IL having smaller

anions.21 We attributed differences in DC between pyr13FSI and pyr13TFSI to the formation of a more compact inner layer for the former system, and a faster displacement from the innermost layer of O(FSI) than of O(TFSI). For the FSI system we showed that the U-shaped character of DC is inherited from the innermost layer and it is caused by a strong adsorption of O(FSI) to the surface near PZC. We report a surprisingly large electrode potential window where the DC decreases steeper than the asymptotic limit U-0.5 predicted for simple EDL models. Simulations at large potentials, i.e., up to 50 V for this system, revealed that when the inner layer becomes overcrowded and grows in thickness, the DC clearly obeys the U-0.5 asymptotic dependence predicted by Kornyshev.9 Finally, we note that, although large induced electrostatic fields can relax in picoseconds via fast vibration modes,48,49 the electrode discharge upon a sudden change in the potential and the EDL formation occurs over time scales of hundreds of picoseconds.

’ ASSOCIATED CONTENT

bS

Supporting Information. Development and validation of the united atom force field, several additional profiles within double layers, and several snapshots of the electrolyte arrangement near electrode. This material is available free of charge via the Internet at http://pubs.acs.org.

’ ACKNOWLEDGMENT The authors are grateful to the Department of Energy for contract grant DE-SC0001912 to the University of Utah. The united atom force field development was supported by the Air Force Office of Scientific Research, Department of the Air Force contract no. FA9550-09-C-0110 to Wasatch Molecular Inc. and University of Utah. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. ’ REFERENCES (1) (a) Seki, S.; Kobayashi, Y.; Miyashiro, H.; Ohno, Y.; Usami, A.; Mita, Y.; Kihira, N.; Watanabe, M.; Terada, N. J. Phys. Chem. C 2009, 113, 6596. (b) Castriota, M.; Caruso, T.; Agostino, R. G.; Cazzanelli, E.; Henderson, W. A.; Passerini, S. J. Phys. Chem. A 2005, 109, 92. (c) Hayashi, K.; Nemoto, Y.; Akuto, K.; Sakurai, Y. J. Power Sources 2005, 146, 689. (d) Garcia, B.; Lavalle’e, S.; Perron, G.; Michot, C.; Armand, M. Electrochim. Acta 2004, 49, 4583. (e) Sato, T.; Maruo, T.; Marukane, S.; Takagi, K. J. Power Sources 2004, 138, 253. 3083


The Journal of Physical Chemistry B (2) (a) Arico, A. S.; Bruce, P.; Scrosati, B.; Tarascon, J.-M.; Schalkwijk W. V. Nat. Mater. 2005, 4, 366. (b) Chmiola, J.; Yushin, G.; Gogotsi, Y.; Portet, C.; Simon, P.; Taberna, P. L. Science 2006, 313, 1760. (c) Tao, F.; Bernasek, S. L. J. Am. Chem. Soc. 2005, 127, 12750. (d) Su, Y.-Z.; Fu, Y.-C.; Yan, J.-W.; Chen, Z.-B.; Mao, B.-W. Angew. Chem., Int. Ed. 2009, 48, 5148. (3) Wang, P.; Zakeeruddin, S. M.; Compte, P.; Exnar, I.; Gratzel, M. J. Am. Chem. Soc. 2003, 125, 1166. (4) Souza, R. F.; Padilha, J. C.; Goncalves, R. S.; Dupont, J. Electrochem. Commun. 2003, 5, 728. (5) (a) Esperanca, M. S. S. J.; Lopes, J. N. C.; Tariq, M.; Santos, L. M. N. B. F.; Magnee, J. W.; Rebelo, P. N. J. Chem. Eng. Data 2010, 55, 3. (b) Apricio, S.; Atilhan, M.; Karadas, F. Ind. Eng. Chem. Res. 2010, 49, 9580. (6) Huddleston, J. G.; Visser, A. E.; Reichert, W. M.; Willauer, H. D.; Broker, G. A.; Rogers, R. D. Green Chem. 2001, 3, 156. (7) Rodgers, E. I.; Sljukic, B.; Hardacre, C.; Compton, R. G. J. Chem. Eng. Data 2009, 54, 2049. (8) (a) Largeot, C.; Portet, C.; Chmiola, J.; Taberna, P. L.; Gogotsi, Y.; Simon, P. J. Am. Chem. Soc. 2008, 130, 2730. (b) Simon, P.; Gogotsi, I. Nat. Mater. 2008, 7, 845. (c) Lust, E.; Kallip, S.; Moller, S.; Janes, A.; Sammelselg, V.; Miidla, P.; Vaartnou, M.; Lust, K. J. Electrochem. Soc. 2010, 150, E175. (9) Kornyshev, A. A. J. Phys. Chem. B 2007, 111, 5545. (10) For example, the Gouy-Chapman-Stern model can predict a minimum in DC at PZC, but it cannot explain the bell-shaped DC experimentally observed for several ionic liquids. (11) Kobrak, M. ECS Trans. 2010, 33, 411. (12) Lamperski, S.; Outhwaite, C. W.; Bhuiyan, L. B. J. Phys. Chem. B 2009, 113, 8925. (13) Alam, M. T.; Islam, M. M.; Okijima, T.; Ohsaka, T. J. Phys. Chem. C 2009, 113, 6569. (14) Su, Y. Z.; Fu, Y. C.; Yan, J. W.; Chen, Z. B.; Mao, B. W. Angew Chem. Int. Ed. 2009, 48, 5148. (15) (a) Langkau, T.; Baltruschat, H. Electrochim. Acta 2002, 47, 1595. (b) Doubova, L. M.; Trasatti, S. Condens. Matter Phys. 2001, 4, 759. (c) Costa, R.; Figueirodo, M.; Pereira, C. M.; Silva, F. Electrochim. Acta 2010, 55, 8916. (16) Islam, M. M.; Alam, M. T.; Okijima, T.; Oshaka, T. J. Phys. Chem. C 2009, 113, 3386. (17) Lockett, V.; Sedev, R.; Ralston, J.; Horne, M.; Rodopoulos, T. J. Phys. Chem. C 2008, 112, 7486. (18) Islam, M. M.; Alam, M. T.; Okijima, T.; Oshaka, T. Electrochem. Commun. 2007, 9, 2370. (19) Figueiredo, M.; Gomes, C.; Costa, R.; Martis, A.; Pereora, C. M.; Silva, F. Electrochim. Acta 2009, 54, 2630. (20) (a) Silva, F.; Gomes, C.; Figueiredo, M.; Renata, C.; Martins, A.; Pereira, C. M. J. Electroanal. Chem. 2008, 662, 153. (b) Costa, R.; Pereira, C. M.; Silva, F. Phys. Chem. Chem. Phys. 2010, 12, 11125. (21) Lockett, V.; Horne, M.; Sedev, R.; Rodopoulos, T.; Ralston, J. Phys. Chem. Chem. Phys. 2010, 12, 12499–12512. (22) Caprio, D.; Stafiej, J.; Borkowska, Z. J. Electroanal. Chem. 2005, 582, 41. (23) Fedorov, M. V.; Georgi, N.; Kornyshev, A. A. Electrochem. Commun. 2010, 12, 296. (24) Georgi, N.; Kornyshev, A. A.; Fedorov, M. V. J. Electroanal. Chem. 2010, 649, 261. (25) Vatamanu, J.; Borodin, O.; Smith, G. D. J. Am. Chem. Soc. 2010, 132, 14825. (26) Trulsson, M.; Algotsson, J.; Forsman, J.; Woodward, C. E. J. Phys. Chem. Lett. 2010, 1, 1191. (27) (a) Zarzycki, P.; M. Rosso, K., M. J. Phys. Chem. C 2010, 114, 10019. (b) Zarzycki, P.; Kerisit, S.; Rosso, K. M. J. Phys. Chem. C 2010, 114, 8905. (28) (a) Yang, L.; Fishbine, B. H.; Migliori, A.; Pratt, L. R. J. Am. Chem. Soc. 2009, 131, 12373. (b) Wander, M. C. F.; Shufor, K. L. J. Phys. Chem. C 2010, 114, 20539. (29) Tazi, S.; Salanne, M.; Simon, C.; Turq, P.; Pounds, M.; Madden, P. A. J. Phys. Chem. B 2010, 114, 8453.

ARTICLE

(30) (a) Yang, L.; Fishbine, B. H.; Migliori, A.; Pratt, L. R. J. Am. Chem. Soc. 2009, 131, 12373. (b) Shim, Y.; Kim, H. J. ACS Nano 2009, 3, 1693. (c) Shim, Y.; Kim, H. J. ACS Nano 2010, 4, 2345. (c) Singh, R.; Monk, J.; Hung, F. R. J. Phys. Chem. C 2010, 114, 15478. (d) Shi, W.; Sorescu, D. C. J. Phys. Chem. B 2010, 114, 15029. (31) Borodin, O. J. Phys. Chem. B 2009, 113, 11463–11478. (32) Borodin, O.; Gorecki, W.; Smith, G. D.; Armand, M. J. Phys. Chem. B 2010, 114, 6786–6798. (33) Bedrov, D.; Borodin, O.; Li, Z.; Smith, G. D. J. Phys. Chem. B 2010, 114 (15), 4984–4997. (34) Reed, S. K.; Lanning, O. J.; Madden, P. A. J. Chem. Phys. 2007, 126, 084704. (35) Vatamanu, J.; Borodin, O.; Smith., G. D. Phys. Chem. Chem. Phys. 2010, 12, 183. (36) Gingrich, T. R.; Wilson, M. Phys. Chem. Phys. 2010, 500, 178. (37) Siepmann, J. I.; Sprik, M. J. Chem. Phys. 1995, 102, 511. (38) (a) Kawata, M.; Mikami, M. Chem. Phys. Lett. 2001, 340, 157. (b) Kawata, M.; Nagashima, U. Chem. Phys. Lett. 2001, 340, 165. (c) Kawata, M.; Mikami, M.; Nagashima, U. J. Chem. Phys. 2002, 116, 3430. (d) Kawata, M.; Mikami, M.; Nagashima, U. J. Chem. Phys. 2002, 117, 3526. (d) Kawata, M.; Mikami, M.; Nagashima, U. J. Chem. Phys. 2001, 115, 4457. (39) Tuckerman, X. M. E.; Berne, B. J.; Martyna., G. J. J. Chem. Phys. 1990, 97, 1990. (40) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (41) Nanbu, N.; Sasaki, Y.; Kitamura, F. Electrochem. Commun. 2003, 5, 383. (42) Nanbu, N.; Kato, T.; Sasaki, Y.; Kitamura, F. Electrochemistry 2005, 73, 610. (43) (a) Feng, G.; Huang, J.; Sumpter, B. G.; Meunier, V.; Qiao, R. Phys. Chem. Chem. Phys. 2010, 12, 5468. (b) Feng, G.; Zhang, J. S.; Qiao, R. J. Phys. Chem. C 2009, 113, 4549. (c) Wang, S.; Li, S.; Cao, Z.; Yan, T. J. Phys. Chem. C 2010, 114, 990. (d) Feng, G.; Qiao, R.; Huang, J.; Sumpter, B. G.; Meunier, V. ACS Nano 2010, 4, 2382. (44) Fedorov, M. V.; Kornyshev, A. A. J. Phys. Chem. B 2008, 112, 11868. (45) Fedorov, M. V.; Kornyshev, A. A. Electrochim. Acta 2008, 53, 6835. (46) Lauw, Y.; Horne, M. D.; Rodopoulos, T.; Leermarks, F. A. M. Phys. Rev. Lett. 2009, 103, 117801. (47) Bazant, M. Z.; Storey, D. B.; Kornishev, A. A. arXiv:1010.3490v2, 2011. (48) Pinnila, C.; Del Popolo, M. G.; Kohanoff, J.; Lynden-Bell, R. M. J. Phys. Chem. B 2007, 111, 4877. (49) Lanning, O.; Madden, P. J. Chem. Phys. 2004, 108, 11069. (50) Urahata, S. M.; Ribeiro, M. C. C. J. Chem. Phys. 2006, 124, 074513. (51) (a) Bhargava, B. L.; Balasubramian, S. J. Chem. Phys. 2005, 123, 144505. (b) Triolo, A.; Russina, O.; Hardacre, C.; Nieuwenhuyzen, M.; Gonzalez, M. A.; Grimm, H. J. Phys. Chem. B 2005, 109, 22061. (c) Del Popolo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 1744. (d) Margulis, C. J.; Stern, H. A.; Berne, B. J. J. Phys. Chem. B 2002, 106, 12017.

3084


Molecular Simulations of the Electric Double Layer Structure

Recommend Documents