Capacitance Optimization in Nanoscale ... - ACS Publications

Jul 9, 2015 - supercapacitors comprising graphene slit electrodes in ... pores because of the interference between internal double layers (largely the...
0 downloads 0 Views 7MB Size
Download PDF
Article pubs.acs.org/JPCC

Capacitance Optimization in Nanoscale Electrochemical Supercapacitors Srinivasa Rao Varanasi and Suresh K. Bhatia* School of Chemical Engineering, The University of Queensland, Brisbane, QLD 4072, Australia S Supporting Information *

ABSTRACT: We perform constant voltage Gibbs ensemble based grand canonical Monte Carlo simulations for nanosized supercapacitors comprising graphene slit electrodes in symmetric and asymmetric electrolytes. Our simulations demonstrate that external electrolyte at the electrode surface can be exploited to positively influence the structure and packing of that inside the slit, when the system is engineered to allow these to interact. Oscillatory dependence of capacitance on slit-pore size, seen in recent results from molecular dynamics simulation and density functional theory, is observed also in our Monte Carlo simulations. A detailed analysis suggests that maximum in capacitance occurs in subnanometre pores because of the interference between internal double layers (largely the Helmholtz parts) on the opposite sides of the slit, expelling the co-ions; and that the oscillatory character of capacitance with pore width is due to relative changes in counterion and co-ion populations with pore width, also dictated by the interference process between the two internal double layers. Our simulations with size-asymmetric and size-symmetric electrolytes with different sets of electrode pairs reveal that when the pore widths of both the electrodes are close to their respective counterion sizes, the electrodes store maximum charge density, yielding maximum capacitance. Thus, it is demonstrated that for asymmetric electrolytes optimum capacitance is obtained using a correspondingly asymmetric electrode combination.

1. INTRODUCTION Electrochemical double layer supercapacitors (EDLCs) are an important class of energy storage systems that exploit the high surface area of porous carbon structures in order to store energy within electric double layers. Such supercapacitors have received considerable attention because they possess high power density and long life, while being environmentally friendly.1 Besides high surface area, efficient EDLC supercapacitors require materials with pores that are adapted to the size of the electrolyte ions. Both experimental and theoretical investigations have revealed anomalous behavior of capacitance with respect to nanopore size and applied voltage; this has been the subject of much attention with respect to optimizing supercapacitor materials for improved performance.2−7 Experimental measurements of Chmiola et al. and Largeotet al.,8−11 using Ti−CDC electrodes and ethyl methylimidazolium bis(trifluoromethanesulfonyl) imide (EMI-TFSI) ionic liquid as well as solvated ions, indicate that the specific capacitance (per unit surface area) exhibits counterintuitive behavior with variation in average pore diameter8−11the capacitance increases when the pore size falls below 1 nm, reaches a maximum at a pore size equal to the characteristic size of the ions, and decreases sharply at even smaller pore sizes. In contrast, it has also been found that the specific capacitance is independent of the pore size,12 based on studies with several porous carbons with varying structures in quaternary © XXXX American Chemical Society

ammonium electrolytes (TEA−BF4/acetonitrile). In explaining this difference, it has been observed by Centeno et al.12 that the surface areas used by Chmiola et al.8 correspond to the BET area, which may be inaccurate and not representative for micropores, at least for pores smaller than 1 nm. However, Chmiola et al.8 have found the same counterintuitive behavior in capacitance even with the surface area estimated by the NLDFT (non-local density functional theory) approach (Supporting Information of ref 8). A more likely explanation is that the carbons used by Centeno et al.12 have wider pore size distributions covering a larger range of pore sizes, as several porous structures were involved in the study, potentially smearing out any oscillatory behavior with pore size. Therefore, the precise reason for the discrepancy between these two results is not immediately clear, and such conflicting results have stimulated much theoretical activity in order to gain deeper understanding of the underlying mechanisms. Theoretically, molecular dynamics (MD) and Monte Carlo (MC) simulations are now established as powerful techniques to investigate system behavior at the nanoscale, as they consider the intermolecular interactions that dictate system response; such details are generally absent in continuum models. A Received: May 4, 2015 Revised: June 15, 2015

A

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

counterion, in order to obtain maximum capacitance. We have performed simulations on both size asymmetric as well as size symmetric electrolytes in combination with different electrode pairs in order to determine the optimum pore sizes for a given electrolyte and vice versa to obtain maximum capacitance.

common observation from both MD and MC studies in the literature is that of increase in capacitance with decrease in pore width for subnanometer pores.13−16 In addition, unexpected oscillation of capacitance with pore width has been observed in the case of graphene like planar electrodes,17−22 in MD as well as density functional theory studies The anomalous peak in the subnanometer range is attributed to desolvation of ions in subnanometer pores,15 but this mechanism is not applicable to similar behavior observed in the case of pure ionic liquids.10,19 According to Jianget al.,19,22 and Feng and Cummings,16 the anomalous behavior of capacitance with pore size is a result of overlap between the two opposing double layers in a pore, and the solvent plays a critical role in determining this behavior in the case of solvent based electrolytes.22 Nevertheless, while the concept of “overlapping double layers” appears promising, a detailed analysis is necessary to establish this. In most simulations, the electrodes are modeled as thin infinite slabs placed few nanometers apart, with a bulk region in between. In such a setup, the number of ions at the electrode− electrolyte interface is not negligible as compared to that in the bulk. However, in practice, in a typical supercapacitor, the number of such interfacial ions is relatively small compared to the bulk. It is therefore necessary to use larger system sizes in order to mimic real supercapacitors, with greater separation between electrodes, so that a significant bulk region is available between their external double layers. Such an approach is likely to be computationally prohibitive because of the much larger number of ions in the simulation cell, with considerable part of the effort being expended in simulating the bulk region. An alternative approach is to perform simulations in the grand canonical ensemble in which the system can exchange particles with an infinite external reservoir. While there are a number of reports in the literature based on grand canonical ensemble Monte Carlo simulations,23−26 in all of these studies, as well as others based on MD smulations,15,19,27,28 electrodes are kept few nanometers apart and a large amount of simulation time has been spent in simulating the bulk region in between electrodes, limiting the ability to simulate larger system sizes for longer simulation times. This issue has been recently addressed by Punnathanam,29 who has proposed a grand canonical Monte Carlo method based on the Gibbs ensemble, modeling the electrodes as two independent phases. Nevertheless, as in most prior simulation studies, particularly those considering planar graphene slit-pore electrodes, the external double layers have been ignored, which alters the fluid structure within the electrodes as well as the charge density on the electrodes, and therefore the capacitance. While this effect will be small in macroscopic systems, and comprise an end effect, it will be significant in nanoscale devices and sensors,30,31 in on-chip supercapacitors,32 and in nanoelectromechanical systems (NEMS), an area of growing importance.33,34 Using constant voltage Gibbs ensemble grand canonical Monte Carlo simulations, we demonstrate here the importance of considering the fluid external to the slit pore electrode, and its influence on the properties of the system, in particular the capacitance, for nanoscale devices where end effects have significant influence. We also investigate the anomalous dependence of capacitance on the pore width, and analyze system properties in order to obtain microscopic insight into this anomaly. The fact that electrolytes (solvent-based or ionic liquid) are not always size symmetric (having cation and anions of equal sizes) suggests the need to optimize the pore size of each of electrode to be commensurate with its respective

2. SYSTEM, MODEL, AND METHODS 2.1. Simulation System. Our system comprises two graphene slit pore electrodes, with each slit comprising two parallel graphene sheets separated by a distance equal to the slit width (center-to-center distance), and an ionic liquid electrolyte added to the slit during the grand canonical Monte Carlo simulations and which fills the pore. Each graphene sheet has dimensions of 3 nm × 3.3 nm, and consists of 400 carbon atoms. We have considered external fluid layers outside the slits, illustrated in Figure 1, of thickness (Lze) 0, 1, 2, 3, 4, 5, 6,

Figure 1. Schematic of electrode configuration in Gibb’s ensemble based grand canonical Monte Carlo simulations at constant voltage. Here W is the pore width and Lze is the maximum thickness of the external fluid region on the outer surface of the electrode plates.

and 7 nm, to study the sensitivity of the ion density distributions and the capacitance to the electrode spacing. This allows determination of the minimum thickness that must be considered to capture the external double layer and ensure convergence of the results. Periodic boundary conditions have been imposed only in two directions (x and y). The separation between the electrodes is such that the distance between the edges of their external fluid layers is 5 nm along the nonperiodic direction (z-axis), i.e. the spacing between the electrodes is 2Lze + 5 nm. A constant voltage (potential difference, chosen as 500 mV in this study) is applied to the electrodes, and the charge on each electrode is distributed in the form of point charges located at the centers of the carbon atoms. The intramolecular degrees of freedom in the graphene sheets have been frozen in all the simulations reported here; i.e., the graphene sheets are considered rigid. 2.2. Intermolecular Potentials. All interactions in our system have been treated as a combination of Lennard-Jones (LJ) interactions and a Coulomb term for electrostatic interactions, as follows:

Φtotal =

∑

Φij(rij) (1)

i ,j,i≠j

where ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ qiqj σij σij Φ(rij) = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ + r ⎝ rij ⎠ ⎦ 4πεoεr rij ⎣⎝ ij ⎠

(2)

Here rij is distance between two atoms (or ions) i, j = 1, 2, qi is partial charges on atom (or ion) i, σij, εij are LJ parameters for the i-j interaction, εo is permittivity of vacuum and εr is relative permittivity. Ion−ion interactions in the electrolyte have been modeled considering charged spheres interacting via a combination of B

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Lennard-Jones and Coulomb interactions as described by eq 1. Two different electrolytes have been considered here, namely size symmetric and size asymmetric electrolytes. For the size symmetric ionic liquid, LJ parameters are the same for both cations and anions, having values εii = 0.99 kJ/mol, σii = 0.506 nm and for the size asymmetric ionic liquid, the LJ parameters are εii = 0.99 kJ/mol (both cation and anion), σii = 0.68 nm (cation) and σii = 0.506 nm (anion). The charges on the cation and anion are 0.78e and −0.78e respectively. Charge−charge interactions on the graphene electrodes have been treated considering interaction between point-charges placed at the centers of the carbon atoms comprising the electrodes. LJ interactions between the atoms of the electrodes have not been considered as the intramolecular degrees of freedom have been frozen, and the graphene sheets considered rigid. There are only electrostatic interactions between the electrode charges. Electrochemical dilatometry studies on carbon electrodes have shown that there can be charge-induced microscopic structural changes in the electrode due to pumping of electrons in and out and also due to insertion of electrolyte ions inside smaller pores.35−38 Nevertheless, the slit pore made by rigid graphene sheets in the present study represents an ideal supercapacitor electrode. We have not treated image interactions explicitly in our simulations. Electrode−ion interactions have been treated considering LJ interactions between the electrode carbon atoms and the electrolyte ions in the pore, and electrostatic interactions between the ions and electrode charges located at the centers of the carbon atoms. The LJ parameters for carbons on the electrode atoms have been taken as εcc = 0.231 kJ/mol and σcc = 0.34 nm, and the cross interaction parameters between electrode atoms and ions of the electrolyte have been calculated using standard Lorentz−Berthelot mixing rules. The electrostatic interactions have been calculated using the modified Wolf method.39−44 This method has been demonstrated as a promising alternative to Ewald based techniques, successfully applied to various complex condensed matter systems such as SiC,41 aluminum silicates42 etc. In this method, a damping function (complementary error function, erfc(αr)), where α is the damping or convergence parameter) has been incorporated in the second term of eq 2 (1/r term) in order to achieve convergence within the cutoff distance (analogous to the real space term in Ewald sum) and a correction term is added, ensuring the charge neutrality within the cutoff sphere. Complete details of this method can be found elsewhere.39,42−44 The value of the convergence parameter (α) is chosen to be 2 nm −1, which ensures convergence of electrostatic interaction energies to those obtained by Ewald summation. All the interactions have been calculated with a cutoff length of 1.5 nm, and all the simulations have been performed with a background having dielectric constant εr = 5, that implicitly mimics a low dielectric organic solvent such as bromobenzene, acetonitrile or chloroform. 2.3. Simulation Procedure. We have performed Gibbs ensemble based constant voltage grand canonical Monte Carlo simulations on nanosized graphene slit pore electrodes for several pore widths, following the recent algorithm of Punnathanam.29 The constant voltage grand canonical partition function is given by ∞

Ξ(μ , Φ, V , T ) =

where N

Q (N , σ , V , T ) =

N

∑ ∑

Q AC(N1A , N1C , N , σ , V , T )

N1A = 0 N1C = 0

(4)

where Q AC(N1A , N1C , N , σ , V , T ) =

(qV )2N N1A ! N1C ! N2A ! N2C!

∫ ..... ∫ e−βU ds2N

(5)

and the probability density distribution is given by f (s 2N , N1A , N1C , N2A , N2C , σ ) ∝

(qV )2N e βμN e βA Φσ e−βU N1A ! N1C ! N2A ! N2C!

(6)

Here μ, Φ, V, and T are chemical potential, applied voltage, volume and temperature, respectively. σ is the charge density on the electrode wall of surface area A, β is (kBT)−1 and q is eβμ. The number of cations and anions in each of the electrodes are N1C, N2C and N1A, N2A respectively, where 1, 2 are the indices of the electrodes. The following condition has been imposed to ensure charge neutrality in the overall system: N1A + N2A = N1C + N2C = N

(7)

There are seven different Monte Carlo moves (formulated using eq 6) performed during the simulation, namely, translation, insertion, deletion, charge plus ion transfer (CPIT), ion transfer (IT), charge transfer (CT), and charge transfer within the same electrode (CTW). In a translation move, a random ion is chosen from a randomly selected electrode and moved by a small distance (determined by the maximum allowed acceptance, 35%). In an insertion (or deletion) move, an electrode is chosen randomly and an ion pair is created at random positions (an ion pair is selected for deletion in a deletion move). In a charge plus ion transfer move (CPIT), a randomly chosen ion (cation or anion) from a randomly selected electrode is transferred to the other electrode, and a charge equivalent to the ion’s charge is added/subtracted to/from the origin/destination electrode, respectively. By doing this, the CPIT move preserves the charge neutrality of each of the electrodes (including the external fluid region). In an ion translation (IT) move, a randomly chosen ion (cation or anion) is transferred to the other electrode, unlike the CPIT move, no charge compensation on the respective electrodes is performed and therefore this move does not preserve the charge neutrality of the individual electrodes. In a charge transfer move (CT), a small amount of charge is transferred between electrodes, distributed nonuniformly between electrode atoms, and in a charge transfer within the same electrode (CTW) move, a small amount of charge is transferred between the plates of a randomly chosen electrode. In each charge move, only 300 of the 800 atoms in an electrode are randomly selected and the charge distributed equally among these. All these moves are accepted or rejected according to Metropolis scheme, and the acceptance probability for any move is

∞

⎧ f ⎫ acc PMCmove = min⎨1, new ⎬ ⎩ fold ⎭

∑ ∑ e β[μN +ΦAσ ]Q (N , σ , V , T ) Aσ =−∞ N = 0

(3) C

⎪

⎪

⎪

⎪

(8) DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C where f new and fold represent probability density distributions (see eq 6) for the states after and before the moves, respectively. Further details of the Monte Carlo moves and their acceptance criterion can be found elsewhere.29 The simulations are initiated starting with empty electrodes. The system is equilibrated for 1000 Monte Carlo (MC) cycles, comprising about 1 million MC moves. The number of Monte Carlo moves in each cycle is equal to the total number of ions in the system, while the minimum number of MC moves in each MC cycle is set to be 100. After equilibration, a production run of about 10 000 MC cycles is performed, with the coordinates of the ions and charges on the electrode atoms stored at the end of each MC cycle for subsequent analysis. All the simulations have been performed at a temperature of 300 K and a bulk activity of 1 nm−3. Bulk simulations were performed at this activity and the bulk density of the electrolyte was found to be 4.68 nm−3, close to, though slightly smaller than, what is expected for an ionic liquid containing imidazolium cations (∼6 nm−3).45 Further, the peak positions and peak heights of the radial distribution function qualitatively matched those of Roy et al.,46 as shown in the Supporting Information. The voltage has been kept constant at 500 mV. All Monte Carlo moves have been given equal priority.

Figure 2. Snapshots of cathode (top) and anode (bottom) for a 0.75 nm slit. The electrodes atoms are shown in cyan color. The electrolyte ions are shown in orange (cations) and blue (anions) color. Lze represents the thickness of the external electrolyte region.

3. RESULTS AND DISCUSSION 3.1. Influence of External Double Layer. In the literature, a number of studies13,47−49 have investigated the capacitance in slit pore electrodes, using single layer graphene sheets as pore walls. However, such studies have overlooked the effect of the external double layer and its interaction with the intrapore ions. While in practice, electrodes will comprise a number of stacked graphene sheets, consideration of electroneutrality will require treatment of the external double layer to establish the charge accumulated in the electrode pore walls. Ordinarily, this will be an end effect that is frequently negligible in macroscopically sized electrodes; however, it will be important in systems comprising a single or very few juxtaposed nanochannels, such as may be expected in nanoscale devices and nanoelectromechanical systems (NEMS).41,42 It is therefore pertinent to investigate this effect, and its influence on the capacitance in a single nanosized pore. In including the external layer, only a finite amount of external electrolyte can be included in the simulation, and the question of how much electrolyte to consider then assumes importance. Figure 2 depicts snapshots of the electrodes (cyan spheres) surrounded by size symmetric electrolyte (cations: orange spheres, anions: blue spheres) for a 0.75 nm slit surrounded by a 6 nm layer of external fluid on the both sides of slit (Lze = 6 nm). A single layer of ions can be seen inside the slit where most of the ions are counterions with a few co-ions. The ions are randomly distributed far away from the electrode plates, as expected for a bulk environment. In order to explore the sensitivity of the system properties to the external fluid outside the electrodes and its thickness, we have performed simulations by varying the value of Lze between 0 and 7 nm. Figure 3 illustrates the density profiles of counterions along the nonperiodic direction (z-axis) of the system for different values of Lze between 0 and 7 nm. Two strong split peaks of counterions close to walls are found inside the slit even when there is no external liquid considered (Lze = 0 nm, Figure 3a). The ionic configuration for Lze = 0 nm, illustrated in Figure 4, indicates that only a single layer (somewhat disordered) of counterions and few co-ions is located inside the slit. When Lze

Figure 3. Density profiles of counterions along z-direction for a 1 nm wide slit-pore electrode, for various thicknesses of external fluid on each side. Insets a−c show magnified views of the peaks inside and on the left and right sides of the slit, respectively.

Figure 4. Ion configurations inside and outside the cathode for various thicknesses of external fluid. Legends and color schemes are as indicated in Figure 2.

is increased to 1 nm, the height of the peaks within the slit electrode increases from 11 to 14 nm−3, and two small peaks of D

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

density external electrolyte layers (diffusive and/or bulk layers). The co-ion density in the slit is nearly the same regardless of the value of Lze. From the above observations, we see that both structure and number density of the intrapore electrolyte are sensitive to Lze (i.e the amount of external electrolyte). It may be expected that this effect of Lze on ion density distribution will influence charge density as well as capacitance. To investigate this, differential capacitance has been computed from the simulation results, through fluctuations in charge density, following:29,50,51

counterions are also found at the electrode boundaries (within the dashed circles in Figure 3). Visual inspection of the snapshot corresponding to Lze = 1 nm in Figure 4 also indicates that there are single layers of counterions on either side of the slit and a broader layer of counterions and co-ions, consistent with the two peaks found Figure 3. These two peaks belong to the respective double layers (DLs) outside the slit (Helmholtz and diffusive layers). When Lze is increased to 2 nm, the external peaks at around ±1.5 nm are reduced and a uniform distribution can be seen beyond the peaks, resembling a bulklike environment. When Lze is further increased, from 2 to 7 nm (curves for Lze= 3 nm, 4 and 5 nm are not shown in Figure 3, for clarity), there is a negligible change in the heights of the peaks both inside (Figure 3a) and outside the slit (Figure 3, parts b and c). Thus, convergence of the system and the external layer is achieved by about 3 nm thickness of the external fluid, and only additional fluid is being considered without significant changes in the density of the electrolyte in the system upon increasing Lze from 3 to 7 nm. The snapshots corresponding to Lze = 2 and 4 nm in Figure 4 indicate that only a homogeneous mixture of counterions and co-ions (bulk) is included upon increasing Lze from 2 to 4 nm. Figure 5 depicts the variation of number density of counterions, and co-ions in the slit pore, and the overall

Cd =

⎛ ⟨σ 2⟩ − ⟨σ ⟩2 ⎞ ⎛ ∂σ ⎞ ⎜ ⎟ = A⎜ ⎟ ⎝ ∂V ⎠ kBT ⎠ ⎝

(9)

Here Cd is the differential capacitance, A is the total surface area available to the electrolyte, kB is Boltzmann constant, T is temperature, σ is charge density, and < > represents ensemble average. The charge density σ was found to follow a Gaussian distribution, as shown in the Supporting Information. Figure 6

Figure 6. Variation of capacitance with thickness of the external fluid region, Lze.

depicts the variation of area normalized capacitance as a function of Lze. There is a steep increase in capacitance from 0.075 to 0.127 F/m2 when Lze is increased from 0 to 1 nm. The capacitance increases to 0.134 F/m2 upon inclusion of another 1 nm (Lze = 2 nm) of external fluid. This steep change is due to addition of external double layers which mostly comprise counterions (Helmholtz layers) and a few co-ions. On further increase in Lze up to 7 nm, the capacitance remains nearly constant at 0.13 F/m2, and essentially only additional bulk layers are included in this increase. From the above observations, it is clear that appropriate amount of electrolyte outside the electrode must be considered in order to obtain the correct properties of the supercapacitor system. Furthermore, the external fluid has a positive influence on the capacitance, leading to a doubling of that obtained in its absence (i.e., for Lze = 0 nm). 3.2. Effect of Pore Width on Capacitance. Having established the minimum size of the external fluid region necessary to ensure convergence of the results, simulations were conducted to investigate the pore size dependence of capacitance. Initially we performed simulations using a size symmetric electrolyte with electrodes of different pore widths between 0.70 and 2.1 nm (center to center distance between opposing pore walls), while including 6 nm of fluid outside the slit electrodes (Lze = 6 nm). The pore sizes of cathode and anode are taken to be equal (symmetric electrode pairs), and an applied voltage of 0.5 V is assumed. Figure 7 shows the

Figure 5. Thickness (Lze) variation of number density of counterions (filled black circles) and co-ions (filled red squares) in the slit, and total number density (filled blue rhombuses).

number density of the system, with thickness of the external fluid region, Lze. For Lze= 0, the density of counterions is about 6.5 nm−3 and that of co-ions is close to 4.5 nm−3 and the overall number density is just the sum of these two contributions (11 nm−3). When Lze is increased to 1 nm, the overall number density drops below the counterion density inside the slit due to addition of external electric double layers (Helmholtz and diffusive layers), highlighted by the encircled regions in Figure 3. The overall number density continues its fall up to Lze = 4 nm, and does not vary further as the external liquid region increases in thickness to 7 nm. The counterion density inside the slit remains nearly constant between 2 and 7 nm, and therefore the fall in the overall density is due to addition of low E

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

characteristics in this figure. (i) There are two distinct maxima in capacitance located at 0.75 and 1.5 nm, (ii) there is a steep increase in capacitance when the pore size falls below 1 nm, (iii) capacitance decreases sharply when pore size is decreased below 0.75 nm, (iv) capacitance increases when pore size is increased from 1 to 1.5 nm, and (v) a gradual decrease in capacitance is observed for pore sizes above 1.5 nm. The values of capacitance from our simulations are somewhat higher than those reported in the literature,10,16 most likely because of the inclusion of the external double layers, and because in our simulations we have considered electrolyte ions as single spheres. Another possible reason for the difference is that experiments are done using nanoporous carbons which possess wide pore size distributions rather than a uniform pore width. Nevertheless, the position of capacitance maximum coincides with previously reported results.16 The oscillatory capacitance represents an intriguing class of behavior that has also been observed in molecular dynamics simulations and classical density functional theory calculations.16,19,21 However, while the sharp peak at small pore size has been reported in

Figure 7. Variation of area normalized capacitance as a function of slit width, for applied voltage of 0.5 V.

variation of capacitance as a function of width of the slit pore between 0.7 and 2.1 nm. There are several interesting

Figure 8. Density profiles of ions inside cathode along z-axis for 0.7, 0.75, 1.5, and 2.1 nm (center- to- center distance) slit pores. F

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C experiments,8,10,11 the oscillatory behavior is not realized in practice as the carbons used (TiC-CDC, SiC-CDC) have a distribution of pore sizes, rather than a single pore size. Nevertheless, on the basis of the simulations, the capacitance anomaly cannot be ignored and persists in experiments, despite being calculated from overestimated specific surface area (BET surface area). For deeper insight into the anomalous trend seen in capacitance as a function of pore size, the structure of the ions in the slit pores was investigated. Figure 8 depicts the density profiles of the ions both inside and outside the cathode, for 0.7, 0.75, 1.5, and 2.1 nm slit pores. The density of the ions inside the 0.7 nm pore is almost zero showing that this pore cannot be accessed by electrolyte ions (of size 0.506 nm), evident also from the low value of capacitance. This low capacitance is exclusively contributed by external double layers, similar to an isolated electrode plate immersed in an electrolyte. A single large peak of cations and negligibly small peak of anions can be seen inside the 0.75 nm slit pore cathode, indicating that the 0.75 nm slit can accommodate only a single layer of ions (largely counterions). Two large cation peaks close to the electrode plates and an anion peak at the midplane are found in the 1.5 nm slit. For the 2.1 nm slit, two large cation peaks close to the electrode plates and a small cation peak at the midplane, and two small anion peaks and an anion minimum inside the pore are evident. The external ion density distributions outside the slit are nearly similar for all the slit widths, and resemble those observed for free planar surfaces.52,53 Snapshots of electrolyte configurations in the cathode are depicted in Figure 9. It is seen that only a single counterion

Figure 10 depicts the charge density contributions from both cations and anions inside the slit pore cathode, as a function of

Figure 10. Variation of charge density of counterions and co-ions (cations and anions respectively, in the case of cathode) with pore width.

pore width. The cation contribution increases somewhat more steeply than that of the anions as pore width increases from 1 to 1.5 nm, but on further increase in pore width the cation contribution remains constant, while that of the anions continues to increase. This results in a small peak in the effective charge density at 1.5 nm. A decrease in the separation from 1.5 to 1.0 nm induces a steep reduction in both cation and anion charge densities leading to a minimum in effective charge density at 1.0 nm. Further reduction in the pore width, from 1 to 0.75 nm, leads to a sharp increase in the cation charge density and a steep decrease in anion charge density. The steep increase/decrease in cation/anion charge density is due to attractive/repulsive interaction from the two graphene surfaces respectively, and at 0.70 nm slit width, the contribution from both the ions is zero as neither can enter the pore. The effective charge density at 0.75 nm is largely contributed by the counterions (a single layer of cations, Figure 8), also evident from the single large peak in the density profile (Figure 8). The subsequent decrease in the charge density with increase in pore width (up to 1.0 nm) is due to entry of co-ions in the pore, the increase from 1.0 to 1.5 nm and the decrease from 1.5 to 2.1 nm can be attributed to the relative changes in counterion and co-ion populations. From the above observations, it can be surmised that the decaying oscillatory character of charge density variation with pore width is a result of interference between charge density profiles of two isolated graphene plates (anologous to the external profiles shown in Figure 8) when brought close to each other. If the pore width is such that only the Helmholtz layers (counterion peaks close to the graphene plates), then there will be a constructive interference between the two Helmholtz layers of the two opposite graphene plates, leading to a dominating counterion population. Any increase in pore width will lead to an inclusion of co-ion layers in the inteference process, implying a decrease in effective charge density. As a result, this interference process is governed by the extent to which the individual double layers (of graphene plates) can fit into the slit. The changes in charge density distributions in the slit with pore width will be reflected in the total energy of the system, in particular contribution from interactions between ions and electrode plates.

Figure 9. Snapshots of electrolyte and electrode configurations in cathode, for pore widths (a) 0.70, (b) 0.75, (c) 1.5, and (d) 2.1 nm. The color scheme is as indicated in Figure 2.

(orange sphere) enters the slit when the slit width is 0.7 nm (Figure 9a). A single layer of counterions is found inside the 0.75 nm slit (Figure 9b), in accordance with the single peak shown in Figure 8. In case of the 1.5 nm slit (Figure 9c), there are two counterion layers close to the slit walls and a co-ion layer is found at the center of the slit, representing two sharp cation peaks and a small anion peak, shown in Figure 8. Inside the 2.1 nm slit (Figure 9d), there are two counterion layers close to the slit walls, two co-ion layers next to the counterion layers and a counterion layer in the middle are found. The ionic configuration in the external layers far away from the slit is homogeneous, resembling bulk-like structure. G

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

vice versa, the charge density and therefore the capacitance will have a similar trend as ion charge density against pore width. 3.3. Optimization of Electrode Pore Widths for Maximum Capacitance. From the above results for size symmetric electrolyte with size symmetric electrodes (pore sizes of cathode and anode are equal), it is clear that when the pore size is such that only a single layer of counterions can be accommodated then the capacitance exhibits a maximum. However, real electrolytes (either ionic liquids or solvent based) are not always symmetric as far as sizes of the ions are concerned and therefore it is likely that asymmetric electrodes having pores of size adapted to the corresponding ion will maximize the capacitance. To investigate this, we have considered an asymmetric electrolyte (cation, 0.68 nm; anion, 0.506 nm) and a symmetric electrolyte (0.506, 0.506 nm) with three different combination of pore sizes for the electrodes. Noting that the pore sizes of a given electrode should be commensurate with the size of the corresponding counterion for maximum capacitance, we have considered pore size combinations of (0.75 nm, 0.75 nm), (1.0 nm, 0.75 nm) and (1.0 nm, 1.0 nm) respectively for (cathode, anode). The values of area normalized differential capacitance for both asymmetric and symmetric electrolytes with the above electrode combinations are listed in Table 1. Errors were estimated by dividing the simulation run into 10 blocks and calculating the standard deviation with respect to the average; these are reported in parentheses. For the symmetric electrolyte, the capacitance with the (0.75 nm, 0.75 nm) electrode pair is about 3.5 times higher than that of the other two electrode combinations whereas the capacitance for the asymmetric electrolyte with the (1.0 nm, 0.75 nm) electrode pair is about 1.5 times higher than for the other two electrode pairs. As shown earlier, maximum capacitance occurs in the case of symmetric electrolyte with a symmetric electrode pair (0.75 nm, 0.75 nm). This is due to constructive interference between the two double layers (largely from the Helmholtz layers) of the individual graphene plates, resulting in a dominating population of counterions inside the slit, and leading to high capacitance. The co-ions will also be accommodated inside the slit when the width of any of the electrodes is increased. The presence of co-ions will lead to a decrease in capacitance by neutralizing part of counterion charge. Also, the capacitance is low when only one of the electrode pores is commensurate with its counterion size and the other is larger than its respective counterion. For instance, for the (1.0 nm, 1.0 nm) electrode pair, the cathode is commensurate with the cation, while the space available in the 1.0 nm anode pore is significantly larger than the anion size of 0.506 nm. The maximum capacitance in the case of the size symmetric electrolyte is higher than for the asymmetric case because the volume occupied by the cations of the symmetric electrolyte (0.506 nm) is lower than that occupied by those in

Figure 11 depicts the two energy contributions from ionelectrode interactions for the whole simulation cell plotted

Figure 11. Variation of energy contributed from ion-electrode interactions as a function of pore width. The contributions from Lennard-Jones and Coulomb interactions are shown in filled black spheres and filled red rectangles, respectively.

against pore width. The decrease in Lennard-Jones energy contribution with decrease in pore width below 0.7 nm (black spheres in Figure 11) can be attributed to the reduction in the total number of ions (cations and anions) inside the pore. The slight fall of LJ energy in the subnanometer range (from 1.0 to 0.75 nm) is likely due to packing effects related to the increased number density of counterions. The energy contributions from LJ and Coulomb interactions at 0.70 nm are exclusively from the external ions since there are no ions inside the slit. The Coulomb contributions show a slow decrease as pore width decreases from 2.1 to 1.5 nm, and a slow increase with further decrease in pore width to 1.0 nm. However, in the subnanometer range, the Coulomb contributions show a large decrease in the energy (increase in magnitude) up to 0.75 nm. This change is in accordance with the changes in total ion charge density with pore width. The greater Coulomb energy contribution at 0.75 nm is attributed to the attractive interactions between the counterions and the electrode. It is not screened by co-ion-electrode repulsions as the co-ions are completely expelled from the slit at this pore size. The above results indicate that charge densities of ions that are trapped inside the slit as a function of pore width (Figure 10) supported by the changes in the ion-electrode energy (Figure 11) with pore width. As the ion charge density inside the slit will induce a counter charge (which may not be equal and opposite, as there will be violation of internal charge neutrality at nonzero voltages54−56) on the electrode plates and

Table 1. Area Normalized Differential Capacitance for Asymmetric and Symmetric Electrolytes with Three Different Combinations of Electrode Pore Sizes, at a Voltage of 0.5 Va capacitance (F/m2) electrode combination Wcathode = 0.75 nm; Wanode = 0.75 nm Wcathode = 1.0 nm; Wanode = 0.75 nm Wcathode = 1.0 nm, Wanode = 1.0 nm a

asymmetric electrolyte (σcation = 0.68 nm; σanion = 0.506 nm) symmetric electrolyte (σcation = 0.506 nm; σanion = 0.506 nm) 0.11 (±0.0075) 0.175 (±0.00623) 0.124 (±0.0051)

0.214 (±0.0091) 0.0612 (±0.0082) 0.0576 (±0.0083)

Error estimates are in parentheses. H

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C the asymmetric electrolyte (0.68 nm), permitting higher cation densities in the 0.75 nm cathode. Figure 12 depicts the structure for the size asymmetric electrolyte for all the pore combinations mentioned above. In

Figure 13. Snapshots of electrolyte and electrode configurations of size symmetric electrolyte in the electrode pair (a, b) (0.75 nm, 0.75 nm), (c, d) (1 nm, 0.75 nm) and (e, f) (1 nm, 1 nm). The color scheme is as indicated in Figure 2.

In the case of the (1 nm, 1 nm) pair, there are counterions as well co-ions in each of the electrode slits proportionately (see Figure 13, parts e and f). The presence of co-ions neutralizes the counterion charge and diminishes the capacitance, as evident from the values of capacitance given in Table 1. The asymmetric electrolyte in combination with smaller pores (0.75 nm, 0.75 nm) delivers lower capacitance since the cation is too large (0.68 nm) to be accommodated inside the cathode pore (0.75 nm), even though the anion-anode pair is commensurate. This electrolyte delivers higher capacitance when the pore size of the cathode is increased to 1 nm to be commensurate with the size of the cation (0.68 nm). The capacitance is lower when width of both the electrodes is 1 nm, due to presence of co-ions (cations) in the anode, neutralizing the counterion charge. To provide further insight into the above observations, we have computed the charge density contributed from the ions trapped inside the slit as well as those outside the slit. The effective values of these charge densities for all the systems are given in Table 2. In the case of the symmetric electrolyte (0.506 nm, 0.506 nm) with a combination of (0.75 nm, 0.75 nm) electrode pair, the effective ionic charge density is higher compared to other electrode combinations. The ionic charge density contributed from the ions trapped inside the slit (±0.13 e/nm3) is higher than that from those outside the slit (±0.08 e/nm3), so that the internally adsorbed ions contribute more to the capacitance than those adsorbed externally. The ionic charge density in the cathode drops significantly when the cathode pore size is increased up to 1 nm (with (1 nm, 0.75 nm) electrode pair) due to entry of co-ions inside the pore. The ionic charge density inside the slit is lower in case of (1 nm, 1 nm) electrode pair as compared to that of (0.75 nm, 0.75 nm) pair, also due to presence of co-ions in both the electrodes. The ionic charge density in the cathode is zero when the asymmetric electrolyte is combined with the (0.75 nm, 0.75 nm) electrode pair, because cations cannot enter the cathode pore and therefore the overall capacitance is contributed only

Figure 12. Snapshots of electrolyte and electrode configurations of size asymmetric electrolyte in various electrode pairs. (a, b) (0.75 nm, 0.75 nm), (c, d) (1 nm, 0.75 nm) and (e, f) (1 nm, 1 nm). The color scheme is as indicated in Figure 2.

case of the (0.75 nm, 0.75 nm) electrode pair, there are no ions found in the cathode pore (see Figure 12a), and a single layer of anions is found in the anode (see Figure 12b), indicating that the cation (0.68 nm) cannot enter the pore and the anion is commensurate with the size of the anode slit. A single layer of counterions is seen in each of the electrode pores in the case of the (1 nm, 0.75 nm) pair (see Figures 12, parts c and d), indicating that both the slits are commensurate with their respective counterions, so that maximum capacitance is anticipated. In the case of the (1 nm, 1 nm) pair, there are largely counterions in each of the electrode slits, with fewer coions in the cathode (see Figure 12e) than in the anode (see Figure 12f). This difference in co-ion population can be attributed to the competition between electrostatic and excluded volume interactions. The presence of co-ions diminishes the capacitance, as evident from the values of capacitance in Table 1. The structure of size symmetric electrolyte (0.506 nm, 0.506 nm) in each of the electrode pore combinations are depicted in Figure 13. In the case of the (0.75 nm, 0.75 nm) electrode pair, there are exclusively counterion layers inside each of the electrode slits (see Figure 13, parts a and b), indicating that both the electrodes are commensurate with their respective electrode pores. As discussed earlier, this combination will deliver maximum capacitance (see Table 1). There are mostly counterions and a few co-ions in each of the slits in case of the (1 nm, 0.75 nm) electrode pair. The cathode accommodates more co-ions (see Figure 13c) due to which some of its counterion charge is neutralized. The anode also accommodates few co-ions (cations) by virtue of the global charge neutrality condition imposed on the system (see Figure 13d). I

DOI: 10.1021/acs.jpcc.5b04254 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Table 2. Effective Ionic Charge Density within as Well as Outside Slit Pores, for Asymmetric and Symmetric Electrolytes, for Electrode Pairs of Three Different Pore Size Combinations, at Applied Voltage of 0.5 Va effective ionic charge density (e/nm3) asymmetric electrolyte σcation = 0.68 nm; σanion = 0.506 nm

a

Φ = 0.5 V

cathode

inside slit outside slit total

0.0 (−0.0094) 0.09039 (0.0069) 0.09039 (−0.0025)

inside slit outside slit total

0.09308 (0.0002) 0.0788 (0.0087) 0.17188 (0.0089)

inside slit outside slit total

0.05956 (0.007) 0.05314 (−0.009) 0.1127 (−0.002)

symmetric electrolyte σcation = 0.506 nm; σanion = 0.506 nm

anode

cathode

Pore Width: Wcathode = 0.75 nm; Wanode = 0.75 nm −0.08337 (−0.0083) 0.1308 (−0.0007) −0.00702 (0.0122) 0.0799 (−0.0039) −0.09039 (0.0025) 0.2107 (−0.0046) Pore Width: Wcathode = 1.0 nm; Wanode = 0.75 nm −0.1214 (−0.015) 0.00478 (−0.001) −0.05048 (0.0061) 0.05081 (−0.0009) −0.17188 (−0.0089) 0.055 59 (−0.0019) Pore Width: Wcathode = 1.0 nm; Wanode = 1.0 nm −0.0588 (0.0078) 0.01287 (−0.0035) −0.0539 (−0.0055) 0.03562 (−0.0019) −0.1127 (0.002) 0.0485 (−0.0054)

anode −0.1328 (0.0009) −0.0779 (−0.0055) −0.2107 (−0.0046) −0.01434 (0.001) −0.04125 (0.0009) −0.05559 (0.0019) −0.01427 (−0.0032) −0.03427 (0.0086) −0.0485 (0.0054)

Values of ionic charge densities at 0.0 V are in parentheses.

capacitance exhibits a maximum only when the electrode pore sizes are commensurate with their respective counterion sizes. In particular, for size-asymmetric electrolyte we show that the electrodes must also be asymmetric, with pore size tuned to be commensurate with their specific counterion, for optimum capacitance. In such case, the electrode charge is essentially contributed by the counterions, whose dominant population is due to the constructive interference of the double layers on the two plates of the electrode (largely the Helmholtz part), by expelling the co-ions.

by the ions adsorbed on the outer surfaces of the electrode plates. For the anode, nearly 90% of the charge contribution is from the counterions (anions) adsorbed inside the anode slit, as its pore width (0.75 nm) is commensurate with the size of the anion (0.506 nm) while excluding the co-ion. When the slit width of cathode is increased to 1 nm, the effective ionic charge density increases to 0.093 e/nm3, and is largely contributed by counterions. The total ionic charge density of the electrodes also increased from ±0.09039 to ±0.17188 e/nm3, which is largely contributed by the ions adsorbed in the slit. This increase is clearly reflected in the value of capacitance for this electrode pair (1 nm, 0.75 nm). In the case of the asymmetric electrolyte with (1 nm, 1 nm) electrode pair, the total ionic charge density in the electrodes decreases to 0.1127 e/nm3, due to the entry of co-ions inside the anode slit. The ionic charge densities at zero voltage are also given Table 2 (values in the parentheses). These charge densities are an order of magnitude smaller than those at 0.5 V, essentially due to the much weaker influence of LJ interactions which promote the ion adsorption in the absence of applied bias voltage. In such cases both cations and anions are present in the slit but most of them are paired due to which the charge density is low, as demonstrated in Supporting Information.

■

ASSOCIATED CONTENT

* Supporting Information S

Validity of electrolyte model, calculation of differential capacitance, and pore filling at zero voltage:. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b04254.

■

AUTHOR INFORMATION

Corresponding Author

*(S.K.B.) Telephone: +61 7 3365 4263. E-mail: s.bhatia@uq. edu.au. Notes

The authors declare no competing financial interest.

■

4. CONCLUSIONS Our results using Gibbs ensemble based grand canonical Monte Carlo simulations indicate that inclusion of external electrolyte layers is necessary to obtain ion density distributions, ion density and capacitance correctly in simulations associated with slit pore electrodes in nanoscale capacitive systems. The appropriate minimum length of external electrolyte layers required is the distance until the bulk part of the electrolyte is reached, and is about 3 nm in the present simulations. The capacitance increases anomalously with decrease in pore size in the sub nanometer range (