The renewable energy (RE) is gaining prominence over fossil fuels giving rise to the necessity of energy storage systems for RE integration. Electrochemical storage systems can be deployed at various scales unlike mechanical storage systems. They can be used in large grid-scale, small scale off-grid applications including electric vehicles and consumer electronics. The Li-ion battery provides the best electrochemical properties but is very expensive for large-scale deployment. Lead acid batteries, however, are cheaper but exhibit lower energy and power density (Parfomak, 2012). Supercapacitors are electrochemical storage systems that store energy in two ways, static double-layer capacitance and electrochemical pseudo-capacitance (Christen and Carlen, 2000). Supercapacitors are becoming more and more attractive because of their high-specific power, fast charge acceptance, and excellent cycle life (Smith et al., 2011). According to Maxwell industries, electromechanical pitch systems that contain supercapacitors provide significant advantages to wind turbine systems (Werkstetter, 2015). These advantages over conventional batteries are long life, predictable aging, wide operating temperature range, and being light weight. The main disadvantages associated with supercapacitors are low-specific energy, high self-discharge (which makes them unattractive for stand by storage purposes), and the cost. However, it is seen that supercapacitor prices have continued to decline more rapidly than batteries (Gogotsi, 2014).

In 2006, Chmiola et al. (2006) showed anomalous increase in capacitance at very small pore sizes (<1 nm). They concluded that traditional belief of continuous decrease in capacitance (Shi, 1996; Pandolfo and Hollenkamp, 2006) with pore size is not correct by showing that the capacitance increases with the decrease in pore size of subnanometer pores. Since then, many researchers have tried to explain this behavior. Later, Chimola et al. also showed that increase in the confinement of the ions leads to enhancement in the capacitance (Largeot et al., 2008). Inclusion of the effect of partial desolvation of ions in the current double-layer theory could lead to a better understanding of the behavior of ions in confined environments. Many researchers are trying to understand fundamental aspects of the electric double-layer capacitance and pseudo-capacitance for enhanced storage properties.

The classical double-layer theories have limitations in explaining the behavior at confined regions of double-layer supercapacitors. Continuum level theories for double-layer formation such as the Gouy–Chapman theory, Stern theory are limited to treatment of dilute or moderate concentrations of electrolytes (Kornyshev, 2007; Perkin et al., 2013). Deviation from these theories is observed for dense ionic systems such as ionic liquids (ILs) which are increasingly being used in supercapacitors. They can be used without a solvent, which results in high concentration of ions inside the porous electrodes. Capacitance for these systems has a maximum at 0 potential rather than a minimum which is predicted by Gouy–Chapman–Stern model for low dilute electrolytes.

The main challenge in modeling EDLCs is to simulate the electrode polarizability, which directly affects ionic organization at its interface with the electrolyte. It may also be required to further calculate electron transfer processes between the electrode and electrolyte (Petersen et al., 2012). Several methods have been trialed in the literature to model electrodes (Siepmann and Sprik, 1995; Straus et al., 1995; Calhoun and Voth, 1996; Crozier et al., 2000; Lanning and Madden, 2004; Fedorov and Kornyshev, 2008; Lamperski and Kłos, 2008; Wang et al., 2014). Initial attempts included assigning the individual electrode atoms with fixed partial charges, which summed up to equal to the total electrode charge (this was know as the constant charge approach) (Lanning and Madden, 2004). Siepmann and Sprik (1995) proposed a more realistic approach (using a metal electrode) by representing the atomic charges as diffuse Gaussian charges. During simulations, the charges were varied in response to electrolyte movement, by minimizing the total coulombic energy under the constraint of constant voltage. Other researchers have further applied the method to explore complicated electrode electrolyte interfaces (Reed et al., 2007; Willard et al., 2009; Merlet et al., 2012). For example, Merlet et al. (2012) were able to simulate a highly porous carbon electrode interfacing with IL electrolyte to explore the role of subnanometer pores on supercapacitance in the system. Although realistic, the method is computationally expensive due to additional computational cost due to energy minimization. Voth et al. (Petersen et al., 2012) have recently proposed an alternate method using method of images, which is computationally less expensive. In this work, we present a simpler method to simulate electrode polarization, which also does not involve energy minimization calculations. In this method, the electrode charges are varied randomly, and the new states are either accepted or rejected according to the Boltzmann criterion. Therefore, there is very little added computational cost.

Molecular simulation techniques (Kiyohara and Asaka, 2007a,b; Yang et al., 2009; Kiyohara et al., 2010; Perkin et al., 2013) and molecular dynamics studies (Merlet et al., 2012) successfully provide insights into development of capacitance at the nanoscale and the workings of an electrical double-layer capacitor. Most of the molecular simulations of EDLCs reported till date consider both the electrodes and the electrolyte region in a single simulation box. Those simulations simulate the electrode–electrolyte interface as well as electrolyte in the bulk. Hence, a larger simulation box is required, which increases the computational load. A simulation technique that focuses only on the electrolyte–electrode interface would be beneficial in speeding-up the simulation. Such an approach is shown by Kiyohara et al. (Kiyohara and Asaka, 2007a,b; Yang et al., 2009; Kiyohara et al., 2010), wherein EDLC was modeled as slit pores in constant voltage grand canonical ensemble, which does not require simulation of electrolyte bulk. Furthermore, Punnathanam (2014) showed that the two electrodes can be simulated in two different boxes using Monte Carlo technique in Gibbs ensemble. This method is computationally efficient as it does not require simulation of the bulk electrolyte, enabling quick equilibration at low system sizes. We therefore extend this technique (Punnathanam, 2014) of molecular simulation of supercapacitors by addition of a novel move of charge variation (CV) on electrode atoms.

We have studied mechanism of EDLC using Monte Carlo simulations in Gibbs ensemble. The study consists of simulating the two electrodes in separate simulation boxes, with the positive electrode in box I and the negative electrode in box II (Figure 1). The Metropolis algorithm consisting of various Monte Carlo moves is derived from the methodology by Punnathanam (2014). The five moves used are translation, insertion of ion pairs, deletion of ion pairs, charge plus ion transfer, and electrode CV; all of them are explained further in this paper. Certain other moves such as simultaneous charge and ion transfer between the two electrodes/boxes are not possible as they violate charge-neutrality requirement in a box. The underlying reason for such a condition is that the two electrodes are essentially charge-neutral at large separations, as shown by Punnathanam (2014).

The partition function (Kiyohara and Asaka, 2007a,b; Limmer et al., 2013) for EDLC in the constant voltage grand canonical ensemble, Ξ, is given by the following equation: Ξμ,Φ,V,T=∑σ=−∞∞∑N=0∞ e−βμN+ΦAσQN,σ,V,T, QNA1,N,σ,V,T=qV2NNA1!NC1!NA2!NC2!∫⋯∫e−βΔUds2N, where NA(1) and NC(1) are the number of anions and cations in box I, respectively. Similarly, NA(2) and NC(2) are in box II. The total number of anions or cations is denoted by N. Hence, NA(1)+NC(1)+NA(2)+NC(2)=N.

Each simulation box is electrically neutral, i.e., NA(1)−NC(1)=NA(2)−NC(2)=Aσ.

The Monte Carlo moves used in the simulations are explained below. Their acceptance criteria are derived by application of detailed balance. An electrode is chosen randomly for each step following Monte Carlo moves.

We modeled the electrode polarizability with Monte Carlo simulations in Gibbs ensemble by implementing the CV move (Punnathanam, 2014) and studied its effect on the electric double layer. The system is simulated at temperature T* = 3.0 and activity Z = 10⁻³ for all the three electrode geometries. The results can be analyzed into three parts as, the effects of confinement/pore size, voltage, and CV moves on the overall capacitance (Figures 3–9).

Geometry (a) is the most confined structure studied. Charge density is expressed as the ratio of reduced charge and the number of electrode atoms and directly affects the capacitance. Figure 3 shows that for this geometry the charge density increased non-linearly with voltage, unlike in other geometries (Figures 5 and 7), wherein the change was linear. The result highlights an important distinction in the light of new understanding of the behavior of highly confined ions in contact with an electrochemical interface. Kornyshev (2007) noted, using a lattice gas model, that, at high packing fraction, the ionic concentration looked more like Fermi distributions rather than Boltzmann distributions, due to the lattice saturation effects. Therefore, the conventional theories such as the Gouy–Chapman theory (Kornyshev, 2007; Perkin et al., 2013) are unable to correctly reproduce the trend at high packing fraction. Furthermore, the electric potential profile very close to electrode interface as predicted by Kornyshev’s model (Kornyshev, 2007), is less steeper than that given by the Gouy–Chapman theory. The analytical model (Kornyshev, 2007) predicted a maximum in capacitance (a bell-shaped curve) at the potential of 0 charge. Therefore, as seen in Figure 4, the model is able to qualitatively reproduce this trend at high ionic packing fraction. The behavior in Figures 3 and 4 is also indicative of the situations when volume effects become important (i.e., effect of confinement). At lower voltages, the ions can readily fill the pore volume, resulting in a higher differential capacitance. However, the differential capacitance decreases with increasing voltage for the geometry as the pore can no longer accommodate additional ions/charges. The packing fraction of ions increases with voltage, resulting in a higher charge density. It is seen that error in estimation of charge density and capacitance decreases with the addition of CV move.

We simulated the system at moderately high packing fraction. This resulted in relatively low acceptance of insertion/deletion moves. An increase in capacitance with electrode pore size (for pores >2 nm in diameter) is usually seen in the literature, due to higher adsorption in larger pores owing to their increased surface area. In our work, the geometries (b) and (c) exhibit this behavior. In these pores, the diffuse layer of the electrolyte ions adjacent to the electrode surface accounted for most of the potential drop. However, in subnanometer pores, in pores whose dimensions are closer to that of the ionic diameters, the ions distort their solvation shells while entering the pores. This effectively reduces the screening due to the solvation shell and increases the net interaction between electrolyte and the electrode, thereby leading to an increased capacitance (Figure 9).

Geometries (b) and (c) have higher pore width (and volume) than geometry (a). It can be seen from the Figures 5–8 that there is no significant effect on charge density and capacitance due to the introduction of the CV move for both these geometries. Hence, it can be concluded that the CV does not influence the capacitance for larger pores or maybe the effect is being nullified by other unknown factors. Further studies are required to estimate these factors.

Monte Carlo simulations of EDLC have been carried out in constant voltage Gibbs ensemble with simple RPM. Values of differential capacitance (with different pore widths) are obtained at different voltages. The decreasing trend of capacitance with increasing voltage in geometry (a) (Figure 4) is attributed to its highly confined nature. A similar trend was also shown by Kornyshev (Kondrat et al., 2011). Figure 4 is also consistent (qualitatively) with Kornyshev’s prediction of bell-shaped variation of capacitance around 0 voltage for ILs (Kornyshev, 2007). Also, the new CV move introduced in the simulation increased the accuracy of calculations. The CV move is simpler to implement and considerably less expensive compared to the iterative energy minimization techniques used in the literature (Siepmann and Sprik, 1995).

The larger pores; geometries (b) and (c) have more space to accommodate additional ions leading to higher capacitance. The electrode structures are basically toy models and hence contribute only qualitatively in the study. However, with the same simulation methodology, one can build more realistic electrode/electrolyte atomic models, by simply incorporating additional Monte Carlo moves for the molecular degrees of freedom. Therefore, the methodology is easily transferable to model complicated systems. Future studies involving such aspects would contribute toward better understanding and design of high performance materials for supercapacitors.

GP, KR and VK have contributed equally to the analysis. PK has helped in drafting the article and critical revision of the work.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.