Computer simulations of ionic liquids at electrochemical interfaces
Ionic liquids are widely used as electrolytes in electrochemical devices. In this context, many experimental and theoretical approaches have been recently developed for characterizing their interface with electrodes. In this perspective article, we review the most recent advances in the field of computer simulations (mainly molecular dynamics). A methodology for simulating electrodes at constant electrical potential is presented. Several types of electrode geometries have been investigated by many groups, in order to model planar, corrugated and porous materials and we summarize the results obtained in terms of the structure of the liquids. This structure governs the quantity of charge which can be stored at the surface of the electrode for a given applied potential, which is the relevant quantity for the highly topical use of ionic liquids in supercapacitors (also known as electrochemical double-layer capacitors). A key feature, which was also shown by atomic force microscopy and surface force apparatus experiments, is the formation of a layered structure for all ionic liquids at the surface of planar electrodes. This organization cannot take place inside nanoporous electrodes, which results in a much better performance for the latter in supercapacitors. The agreement between simulations and electrochemical experiments remains qualitative only though, and we outline future directions which should enhance the predictive power of computer simulations. In the longer term, atomistic simulations will also be applied to the case of electron transfer reactions at the interface, enabling the application to a broader area of problems in electrochemistry, and the few recent works in this field are also commented upon.
Ionic liquids are room-temperature molten salts, composed mostly of organic ions for which almost unlimited structural variations are possible . The scientific and technological importance of ionic liquids derives from a wide-range of applications [2, 3]. For example they can be used as solvents in the textile industry, as lubricants, or as heat-transfer fluids for thermal engines. Due to their excellent properties as solvents they may even enable emerging technologies in such fields. Here we will focus on the electrochemical applications of ionic liquids only. These can be separated in two main families: electro-deposition on the one hand and energy storage and conversion on the other hand. For the former, the large electrochemical windows of many ionic liquids allow for the development of processes that are impossible in water; for example the electroplating of aluminium in order to protect steel from corrosion . As for the very-active field of electrochemical storage of energy, many synthesis routes for novel materials involve the use of ionic liquids . They are also used in replacement of conventional organic solvents as electrolytes in battery , fuel cell or supercapacitor  devices, allowing for the exploitation of their chemical stability across a large electric potential window.
In particular, it is worth emphasizing the case of electrical double layer capacitors (EDLC), which has attracted much attention in recent years . In these devices, the operating voltage is a crucial parameter, which impacts both the power and energy densities . Many groups have therefore tried to develop ionic liquids with large electrochemical windows [11, 12]. The other target properties are a low viscosity, in order to decrease the charge/discharge time, down to very low temperatures , and the capacitance at the ionic liquid/electrode interface. The discovery of nanoporous electrode materials with enhanced capacitive performances when using an ionic liquid with ion sizes matching the pore size, as illustrated in Figure 1, opens the way for a widespread use of supercapacitors in many contexts where high power electrical output is required . As an extension, the behavior of the ionic liquids at charged interfaces even allows one to forecast the development of electroactuators [15, 16] which could be used as artificial muscles, sensors, and even energy generators in turbulent flows or sea-tide.
The surge of experimental results, some of them somewhat unexpected, which has followed this technical interest and been made possible by new experimental methods, has spurred the interest of the theoretical community. Despite the underappreciated early simulation work of Heyes and Clarke , the interface between electrodes and ionic liquids had barely been explored until the last decade. Since then several studies have focused on the structure of ionic liquids on charged flat surfaces [18, 19], and shown that the ion densities exhibit a pronounced oscillatory structure close to the interface. In parallel, the development of a mean-field theory based on the Poisson-Boltzmann lattice-gas model has shown that it is mandatory to account for the finite volume occupied by the ions, resulting in a dramatic departure from the Gouy-Chapman law for the dependence of the capacitance on potential . It is worth noting that the latter work actually preceded experiments, and that numerous studies have reported differential capacitance results in qualitative agreement with the mean-field theory since then [21, 22].
More refined simulation studies have focused on a better understanding of the existence of particular capacitance-potential curve shapes, such as e.g. the “bell-shaped” or “camel-shaped” ones. In particular, the roles of the asymmetry of size between the ions or of the charge distribution inside a given ion have been adressed [23, 24, 25, 26]. It was also shown that for the particular case of the interface between the molten salt LiCl and an aluminium electrode, a potential-induced ordering transition of the adsorbed layer may be at the origin of a sharp peak of the capacitance-potential curve [27, 28]. Such a situation has also been reported for imidazolium ionic liquids systems in experiments based on in situ scanning tunneling microscope [29, 30]. Lately, the adsorption on atomically corrugated electrodes was shown to induce a systematic increase in the differential capacitance compared to the case of flat electrodes [31, 32, 33].
The technologically important case of porous electrode has also recently been investigated. Simulations involving slit-like pores  or carbon nanotubes [34, 35] provided a first insight into the structure of ionic liquids in a confined environment. More recently, the importance of the polarization of the (conducting) electrodes by ions in close proximity to them has been appreciated, and shown to play the key role in allowing the “superionic” state, which accounts for the high degree of charging at the interface inferred from experiments on porous electrodes, to develop. [36, 37]. The mechanism at the origin of the enhanced capacitance in nanoporous electrodes was then fully understood from molecular dynamics simulations involving realistic models of carbon materials and which took into account the polarization in an appropriate way . In order to underline the current vigor of the simulation community in this field, we may mention the following example: the existence of a capacitance varying with an oscillatory behavior depending on the pore size was simultaneously reported by three different groups in the same week [39, 40, 41]!
Computer simulations involving ionic liquids at the molecular scale employ either the molecular dynamics (MD) or Monte-Carlo (MC) method . In both cases, it is necessary to determine the intramolecular and intermolecular interactions at each step of the simulation. Among the various existing methodologies, the size of the systems involved in electrochemical applications limits the options to the use of classical approaches, where all these interactions are defined through an effective force field. The intramolecular and van der Waals (repulsion + dispersion) terms generally differ from one study to another in the analytic expressions involved (e.g. Lennard-Jones vs. Born-Huggins-Meyer for the van der Waals term) or in the parameterization scheme [43, 44, 45, 46], with very little impact on the qualitative results for the structure (the situation is somewhat more complicated for the dynamics, but the latter has barely been addressed in interfacial systems). In some cases, the computational cost is decreased by employing coarse-grained models [47, 48].
More significant are the differences in the treatment of the electrostatic problem and, in particular, the influence of the presence of the charged and conducting electrode on the effective interactions. The electrodes may be treated as ideal conductors and in macroscopic theory the charge induced on the electrode by the applied potential and the charges of the ions in the electrolyte is obtained from the Poisson equation under the condition that, for inside the metal,
where and is the charge density due to the ionic liquid (represented by a set of point charges located on the atomic sites). is the applied potential, which should be uniform inside a conducting electrode. To solve this equation, which must be done “on-the-fly” as the ions move, it is necessary to include specific developments in the simulation codes, and for this reason in most of the simulation studies a much simpler representation of the “electrode” was used. In that case the electrode is either represented by a planar surface with a fixed uniform surface charge density (i.e. with no explicit atomic sites) or with a set of fixed partial charges located on fixed atomic positions. In the remaining text we will designate such simulations as constant charge ones and those in which the induced charge is determined self-consistently by solving equation 1 as constant potential ones.
In parallel, several methods have been developed in order to perform constant potential simulations. In particular, Siepmann and Sprik proposed a model for studying the adsorption of water molecules at a metal surface and on the tip of a model scanning tunneling microscope probe , which was adapted by Reed et al. to the case of electrochemical cells . In this model, the electrode consists of atoms which can be arranged with a crystalline order or in a disordered way; each atom carries a Gaussian charge distribution which has an integrated charge of and is of fixed width :
where is a normalization constant. In our application, the full system consists of two parallel electrodes, which are infinite in the and direction, enclosing a set of melt ions. The Coulomb energy, given by,
has to be expressed through a two-dimensional Ewald summation, for which the correct expressions are provided in reference  (note that the use of three-dimensional Ewald summation may lead to artifacts, although no systematic comparison has been undertaken up to now). By combining equations 1 and 3, we immediately see that the potential experienced by any charge is obtained from the partial derivative of this expression with respect to that charge
The value of the charge on each electrode atom at each time step in a MD procedure is obtained by requiring that the potential experienced by each charge in that electrode be equal to the preset electrode potential value, , thus satisfying the constant potential condition. This condition is achieved by adding a constraint term,
where is the number of electrode atoms, and by minimizing the total potential energy with respect to all the variable charges simultaneously (note that the two electrodes will have different preset potential values and , which will be noted and in the following).
A second approach has also been used in the ionic liquid / electrode context, in which the electrode surface is modeled as an equi-potential smooth surface with a net charge . Here, the entire electrode surface is effectively modeled as a metal foil separating the ions and a set of fixed charges located inside the electrode. The constant electrical potential is enforced on numerical grid points lying over the electrode surface by solving an auxiliary Laplace equation, following a method proposed in another context by Raghunathan and Aluru . This method differs from the Reed et al. approach by treating the surface as smooth and by replacing the self-consistent calculation of the charges by solving Poisson’s equation on a grid.
In their grand canonical MC simulation study of a model ionic liquid in a slit-like metallic nanopore, Kondrat and Kornyshev have employed a third method, which consists in modifying the electrostatic interactions between pairs of ions by introducing an effective screening by the metallic wall . The electrostatic interaction between two ions of respective positions and , charges and , separated by a lateral distance , is given by
where is the dielectric constant of a medium between the plates, is the pore size and is the zero order modified Bessel function of the second kind. The applied electrostatic potential is introduced by adding a shift in the chemical potential of the ions. This approach is limited to pore geometries for which the screening function can be evaluated analytically, and only the first order polarization effects are taken into account, so that the results are mostly qualitative.
More recently, another computationally efficient alternative method was proposed by Petersen et al., in which the polarization of the electrode by the electrolyte is accounted for by introducing explicit image charges together with a fluctuating uniform electrode charge located on the electrode surface. The constant potential condition is then enforced by adding a constant uniform charge . The main advantage is that it avoids the need of a self-consistent calculation at each time step, but it presents the drawback of being a priori limited to the case of electrodes with planar geometries.
It is also worth noting the work by Pádua et al., who have studied the solvation and stabilization of metallic nanoparticles in ionic liquids, an application in which it is not necessary to control the potential [55, 56, 57]. In this work, a polarizable model is used for the nanoparticles, which is done in practice by assigning a freely rotating dipole to each atom.
In practice, constant potential simulations are more easily comparable to experimental data than constant charge ones, because the conditions are more realistic. Furthermore, the first-mentioned method above can be applied to electrode surfaces with arbitrary geometries because the charges are determined self-consistently and can vary from site-to-site depending on the local geometry, whereas in constant charge methods the charges must be preassigned with due regard to the symmetry. In the remainder of the text we will therefore favor the results obtained by using such methods. Nevertheless, constant charge simulations are easier to carry out (no particular modification has to be made in regular simulation packages) and their computational cost is lower (by an order of magnitude). A much larger number of situations have therefore been studied in this framework, and we will also discuss such results when no data from constant potential simulations are available. Note that we have recently shown that, for simple electrode geometries where the average charge distribution can be predicted, in general constant charge simulations provide a fair description of the structure of the liquid at the surface of an electrode, and thus of the capacitance. For the dynamics, the situation is much more dramatic: In particular, upon switching on an applied potential difference, the increase in the temperature, due to the Joule effect, associated with the creation of an electric current across the cell follows Ohm’s law, while unphysically high temperatures rapidly develop when the potential change is induced by changing the charges assigned to each carbon atom.  As a consequence, charging kinetics are largely affected by the simulation method, being unrealistically high for constant charge simulations. [59, 58]
Ii Planar electrodes
The vast majority of the simulation work up to now has been directed toward the understanding of the adsorption of ionic liquids on planar electrodes. In addition to the changes in the simulation method (choice of the force field, of the simulation conditions, of constant charge vs. constant potential, etc), they vary by the nature of the ionic liquids (molten salts [17, 18, 50, 27, 28, 26], imidazolium-based ionic liquids [63, 64, 62, 65], primitive models [23, 66, 24, 25]) and of the electrode (graphite/graphene [63, 62, 65], smooth surface [23, 66, 24, 25], generic metal ,LiFePO  or solid aluminium [27, 28]) involved. Nevertheless the main picture of the structure has remained essentially the same since the pioneering work from Heyes and Clarke : The ionic liquid adopts a layered structure, and the structure of the bulk is recovered after several layers. This particular pattern has since then been confirmed by several experimental techniques such as high-energy X-ray reflectivity , AFM [69, 70, 60] or surface force apparatus (SFA) [71, 61, 72, 73]. Figure 3 illustrates this layering as it is observed in AFM  (top) and SFA  (middle) experiments or in MD simulations  (bottom). Similar patterns are observed for all the ionic liquids: In the illustrated example, the ionic liquids are composed of ions with different sizes and shapes, and the nature of the solid surface also changes, which explains the variations in the distances between successive planes in the different cases. Sum Frequency Generation (SFG) is another experimental technique which is able to probe the vibrational properties of molecules at interfaces. In the case of ionic liquids, Baldelli did not deduce from his experiments the presence of several layers of ions at the surface of an electrode [74, 75, 76]. It is now possible to determine SFG spectra from molecular dynamics simulations with a high degree of accuracy [77, 78], and the application of such methods to ionic liquids interfaces will certainly help to understand this discrepancy in the future.
Within the plane of the adsorbed layers, the ionic liquid can also take a noteworthy organization. Of course, in a given layer, and when the electrode is not charged, the dominant interactions are the short-range repulsion, which controls the packing of ions, and the Coulombic interactions which occur between the various ions. A cation is then surrounded by anions and vice-versa. Nevertheless in situ scanning tunnel microscopy (STM) studies have provided evidence of potential-induced transitions in the adsorbed layers of ionic liquids, with the formation of ordered structures on the surface of the electrode [79, 80, 29, 30]. In computer simulations, Kislenko et al. have first shown that in pure BMI-PF, ions are arranged in a highly defective 2D hexagonal lattice on a neutral graphite electrode, as illustrated from the density contour maps on figure 4. The ordering of the ions in the first adsorbed layers was also observed in the LiCl // aluminium interface , and it was shown that the structure of the ionic liquid film was commensurate with that of the aluminium surface, with different patterning depending on the metallic face oriented towards the liquid . The question of the orientational ordering of the ionic liquid has also been addressed in numerous computer simulations [19, 82, 56, 83, 84, 85, 86].
Upon charging of the electrode, the ionic liquid keeps a layered structure at the interface. The accumulation of charge on the surface of the metal leads to the formation of additional layers: A segregation of the ionic species is progressively observed , which is at the origin of an overscreening effect [18, 23, 66, 87, 88]. The first adsorbed layer carries a charge which is larger in magnitude (but opposite in sign) to that on the electrode,which is then counterbalanced in the subsequent layer. This phenomenon extends in the following layers, until the bulk structure is recovered (i.e. up to several nm). This overscreening effect is maximum for low voltage electrodes, and it disappears when high charging régime is reached [23, 66] (typically for metal surface charge greater than 16.0 C cm, i.e. substantially out of the electrochemical window of any electrolyte, even an ionic liquid !).
The observable property targetted in these simulations is the capacitance of the electrochemical cell. One can define either a differential capacitance, which is determined for each electrode by
where is the average surface charge of the electrode111The symbol is used for quantities which can be calculated for both the positive and negative electrodes and is the potential drop at the electrode/electrolyte interface, or an integral capacitance, given for the full electrochemical cell by
In constant charge simulations, and are imposed and the various potential differences involving , and are obtained from the charge density distribution along the normal to the interface:
where is an arbitrary integration constant (note that this relation is only valid for planar electrodes). In constant potential simulations, and are imposed and equation 9 is also used to determine , but now we have the boundary condition if is inside the left electrode ( in the right electrode). In addition, the surface charges are calculated at each time step as explained in the Methods subsection. Typical fluctuations of the positive electrode surface charge along simulation time are shown on figure 5.
A typical plot, obtained for an electrochemical cell made of graphite electrodes and BMI-PF ionic liquid, is shown on figure 6. In our previous work , we proposed linear fits of the data for both the and regions and observed that this fit was not valid outside the electrochemical stability window. Nevertheless it was not possible to check the existence of a region where the so-called lattice saturation effect occurs, as predicted in the mean field theory developed by Kornyshev . By extending the potential range we can now assess this point; as shown on figure 6, there is a proportional increase of with respect to for potential drops greater than 4 V, which is in agreement with the theoretical prediction. In this régime, there is no more overscreening, and on the contrary the ions cannot pack densely enough inside one layer to counterbalance the electrode charge. Comparison with other simulation work seems to indicate that the potential drop for which lattice saturation effects occur are largely dependent on the nature of the ionic liquid [25, 59, 89], but the differences in the models and methods between the various studies could also be the origin of this variation.
If we now come back to potentials within the electrochemical stability window of the ionic liquids electrolytes, whether there is a general shape for the differential capacitance vs. potential drop (C-V) curves is still an open question. Indeed, many experimental results have been reported following the theoretical work by Kornyshev which predicted the existence of “bell-shaped” curves. Non trivial shapes (e.g. “camel-shape”) were measured [91, 21, 92, 93, 94]. Nevertheless Lockett et al. showed that many hysteresis effects occured in these experiments , an observation which was later confirmed by Drüschler et al. . An important parameter is the structure of the electrode surface, which is often far from the perfect flat surface assumed in the theory since most of the experiments involve either polycrystalline metal surfaces or glassy carbon electrodes . It was also shown that the use of single-frequency measurements could lead to misleading results [98, 99, 100].
Although some of the experimental problems do not arise in computer simulations, where, for example, it is much easier to study a perfect monocrystalline surface than a polycrystalline one, they also have their drawbacks. The most important one is that classical modelling is not suitable for the determination of the potential of zero charge (PZC), which is defined as the value of for which ; of course a value can be extracted from C-V curves but those cannot be compared to experiments because it corresponds to the position of the Fermi level of the electrode, an electronic property which can only be caught at the ab initio level. Ab initio molecular dynamics may provide a solution to this problem in future years , but this technique remains substantially too expensive in computational resources because of the large cells and long relaxation times necessary for the IL interfacial studies – for the moment it is limited to the study of relatively small bulk ionic liquids samples [102, 103, 104]. The surface charges and the relative potentials fluctuate on long timescales so that their averages are associated with substantial error bars; for example on figure 5. In turn, the method used to extract the C-V curve is also subject to numerical uncertainties. For example, the linear fits of the positive and negative branches of the surface charge vs. potential drop plots can be replaced by polynomial fits as suggested by other groups [23, 66, 59]. We are currently exploring a new route to the differential capacitance, namely the equilibrium fluctuations of the electrode charge in constant-potential simulations. The variance of this quantity can indeed be shown to be related to the differential capacitance , which could then be computed without resorting to simulations at different potentials.
Typical C-V plots obtained by Vatamanu et al. from molecular dynamics simulations of the ionic liquid N-methyl-N-propylpyrrolidinium bis(trifluoromethane)sulfonyl imide (pyr-TFSI) near a graphite electrode are shown on figure 7. The authors showed that the various maxima/minima were due to changes in both the ion density and in their orientation towards the surface. The shape of the ions and their asymetry therefore play an important role; for example at 363 K and 393 K the lower value of the maximum of at negative potentials compared to positive potential seems to reflect the relatively high density of the uncharged pyr propyl tails and ring carbon atoms near the graphite surface for the former, while the latter corresponds to an innermost adsorbed layer of highly charged oxygen atoms of the TFSI anion.
Iii Curved/rough electrodes
The study of planar electrodes is important for basic understanding of the electrical double layer formed with ionic liquid electrolytes. But both at the laboratory and industrial scale, EDLC devices employ active matter with a somewhat more complicated structure: for example one can find nanotubes , nano-onions, cluster-assembled nanostructured carbons , carbide-derived nanoporous carbons , etc. In all these setups, the carbon surface is mostly curved, and it is important to understand how this impacts the structure of the adsorbed ionic liquid and the capacitance, inter alia.
Several groups have therefore tackled this problem using molecular dynamics simulations. On the one hand, Feng et al. have simulated, using a constant charge setup, both the outer surface of carbon nanotubes [109, 110] and nano-onions . In both studies, an overall increase of the capacitance per surface area was observed compared to a planar surface. The main result from reference  is shown on figure 8. It reports the variation of the differential capacitance, with respect to the potential drop across each interface (i.e. between the positive electrode and the electrolyte or between the negative electrode and the electrolyte), for a 1-ethyl-3-methylimidazolium bis-(trifluoromethylsylfonyl)imide (EMI-TfN) ionic liquid and three different electrode structures: planar graphite and two nano-onions of different radius (0.705 nm and 1.223 nm). The former shows a peak at 0 V, which does not appear in the two other cases. The lattice saturation occurs at relatively small potentials for the planar electrode (-1.0 V on the cathode side and 0.8 V on the anode side) while it does not occur in the investigated potential range for the nanonions. We therefore expect this effect to play an important role in the observed difference. The authors also showed that when the nano-onion becomes smaller (thus increasing the surface curvature), the number of ions accumulated per unit of area of the electrode rises, which leads to an increase of the capacitance (by 10 % across the whole potential range).
On the other hand, Vatamanu et al. have investigated rough surfaces using constant potential simulations . To this end, they have modelled the carbon electrode by using several graphite layers which are stacked perpendicularly to the interface. They used an ABAB stacking sequence, which leads to a surface with a regular indentation of 1.43 Å due to the relative positions of the atoms at the edges of the A and B layers. The ionic liquid consists of 1-ethyl-3-methylimidazolium cations and bisfluorosulfonylimide anions. In agreement with the work of Feng et al., they showed that a large increase of the capacitance per unit area is obtained, relative to a planar graphene-like surface. The shape of the differential capacitance vs. potential curve is different, though, since a peak is obtained for the roughened interface but not for the planar one. Another important consequence of their study is that even for a simple structure like graphite, the shape of the capacitance vs. potential plot can vary dramatically depending on the surface which is facing the liquid, a point which should be kept in mind when comparing simulations when experimental data. This was also observed in our study of the interface between a LiCl molten salt and two different faces of aluminum electrodes, namely (100) and (110) .
Iv Porous electrodes
As soon as porous electrodes are considered, several additional effects arise. First of all, due to the confinement the ions close to the electrode surface become less coordinated. In recent work, we performed simulations of a supercapacitor consisting of a BMI-PF ionic liquid and carbide-derived carbons (CDCs) electrodes . CDCs are nanoporous carbons with well-defined pore size distribution [108, 111], for which a strong increase of the (integral) capacitance with decreasing pore size was shown experimentally for various electrolytes [112, 8, 113, 114]. In our simulations, we have used structures which were obtained by Palmer et al. using quenched molecular dynamics  and were shown to correspond to CDCs synthesized from crystalline titanium carbide using different chlorination temperatures. Figure 9 displays a snapshot extracted from one of our simulations of this system. For a null voltage difference between the electrodes, the coordination number drops from seven in the bulk to four inside the electrodes .
Then, as soon as a potential is applied between the electrodes, an exchange of ionic species occurs between the electrolyte and those initially absorbed within the nanoporous carbons. In an elegant analytical study, Kondrat and Kornyshev have shown that this is facilitated by the image charges at the surface of the metallic electrode, which exponentially screen the electrostatic interactions between ions close to the surface . Our numerical simulations are in qualitative agreement with this model. The volume occupied by the liquid inside the electrode remains almost constant, but its composition is changed and it is no more electrically neutral ; the net charge at the surface of the electrode exactly balances that of the adsorbed liquid. In addition, the coordination number of anions (resp. cations) inside the negative (resp. positive) electrode decreases even further, passing from 4 to 3, for an applied potential of 1 V. On the contrary, cations (resp. anions) become more highly coordinated inside the pore (their average coordination number is then around 5).
In agreement with experiments, the associated integral capacitance is then very high: In our study, we calculated values of 87 and 125 F per gram of carbon for two different CDC structures and applied potentials of 0.75–1 V, while it was of 30 F g in the case of planar graphite electrodes. In order to understand this difference, it is very useful to compare the corresponding ionic density profiles. Figure 10 displays these profiles for one of the porous carbon structures and for the planar graphite surface. Two key differences are observed. Firstly, the ions approach closer to the surface by 0.7 Å in the porous electrode case. This partly explains the increase of the capacitance, since a shorter distance will give rise to a larger local polarization of the carbon surface. Secondly, due again to the confinement inside the pore, only one layer of ions is now adsorbed on the surface. As a consequence, overscreening cannot take place inside nanoporous carbon electrodes, an effect which is clearly visible when comparing the density of charge on the liquid side of the interface with the surface charge . Although in the case of graphite the charge in the first adsorbed layer reaches a value which is three times higher than the charge of the electrode itself, in the CDC the total charge in the first layer balances exactly that of the electrode; this results in a much better efficiency for the charge storage.
Other recent simulation works have focused on the study of the adsorption of ionic liquids inside nanoporous electrodes. Idealized slit pores or regular carbon structures such as the interior of nanotubes were considered, either by molecular simulations or using a classical density functional theory approach [116, 34, 35, 26, 39, 40, 41, 37, 117, 118]. All these works were in qualitative agreement with the increase of capacitance for pore sizes matching the ionic size observed experimentally and reproduced in figure 1. Interestingly, several studies reported an oscillatory behavior of the capacitance with pore size [39, 40, 41], an effect which could not be observed experimentally because real electrodes used up to now always exhibit a pore size distribution. In addition, our simulations showed that due to the variety of environments encountered by the ions inside the relatively disordered porous carbons, the capacitance is not only linked to the size of the pores but also depends on the local structure inside the porosities , in agreement with the curvature and roughness effects described in the previous section. Finally, it is worth noting that the effect of pore size dispersity has been addressed in a single study so far by Kondrat et al. .
V Concluding remarks
While to date most theoretical work has been devoted to the simulation of pure ionic liquids, most experimental devices actually use ions dissolved in a solvent. Aqueous electrolytes, which are used in batteries and display a high ionic conductivity, have been studied for a long time by experiments, theory and simulations. In fact, the understanding of charged solid/liquid interfaces has been the subject of over a century of theoretical work, from the early studies of Gouy and Chapman and the concept of double layer to the most recent molecular simulations. However aqueous electrolytes offer only a rather narrow electrochemical window and organic solvents such as acetonitrile (ACN) or propylene carbonate (PC) are usually preferred for EDLC applications. Future theoretical work should thus focus on the solvation of organic ions in these liquids and at their interface with electrodes [120, 121]. The challenge will then be to describe how solvation impacts the structure of the electrolyte at the interface and the resulting charge storage process.
Significant theoretical work should also focus on the solid counterpart of the electrolyte, namely the description of the porous electrode. While the importance of using realistic electrode structure to capture the microscopic charging process has been recently demonstrated , resort to simplified theories remains necessary to predict generic trends e.g. of the ratio between molecular and pore sizes or of polydispersity . We thus anticipate that both routes will still provide useful insights. In the former case, the topological analysis of experimental microscopic structures should reveal characteristic properties which influence the capacitance beyond the average pore size or even the pore size distribution. This will improve the understanding of the intimate interplay between the porous structure and the electrolyte, in particular for the dynamical aspect of the charge and discharge. In the latter case, simplified studies will allow for a systematic screening for the design of improved materials.
From the theoretical point of view, it would also be desirable to improve the description of the electrode material. The shift from a constant charge to a constant potential within the classical simulation framework has already been a significant step in the description of the electrode  but it would be useful to have ways of systematically improving the description of the short-range interactions between ions and the atoms of the electrode surface. We note that Pounds et al  implemented a “force-fitting” method to inject first-principles parameterization for a classical ion-electrode force-field of a prescribed functional form - but further systematic investigations are necessary to expose what physical effects should be included in such models, and ab initio simulations of small-scale systems would be useful. We also note that the currently available constant potential methods treat the electrode as an ideal metal. It is not a priori obvious that this is appropriate for a disordered metal  or porous carbon (or whether there are limitations on the appropriate length- or time-scales on which it might be justified), and it is certainly not appropriate for the semiconductor surfaces which one would like to be able to treat.
Several further directions are likely to be explored in the near future within the framework of the currently available methods. The most immediate step, even though not straightforward, is to improve the numerical efficiency of the algorithms to maintain a constant potential condition. Such phenomena as electrowetting at metallic surfaces  or the distribution of ions around particular topological features, like asperities or step-edges are ripe for exploration with the constant potential method. The current review has been confined to non-faradaic process at the electrodes surface, but charge transfer induced by a change in the electrode potential between the electrode and redox species in the electrolyte has been investigated  with the constant potential method within the framework of Marcus Theory. This work has shown how to relate, at the atomistic level, the macroscopic potential applied to the cell to potentials felt by the redox species in the electrolyte.
We are grateful to Alexei Kornyshev, Patrice Simon, Pierre-Louis Taberna, René van Roij and Susan Perkin for useful discussions. We acknowledge the support of the French Agence Nationale de la Recherche (ANR) under grant ANR-2010-BLAN-0933-02 (‘Modeling the Ion Adsorption in Carbon Micropores’). We are grateful for the computing resources on JADE (CINES, French National HPC) obtained through the project x2012096728. This work made use of the facilities of HECToR, the UK’s national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc and NAG Ltd, and funded by the Office of Science and Technology through EPSRC’s High End Computing Programme.
-  J. P. Hallett and T. Welton. 2. Chem. Rev., 111(5):3508–3576, 2011.
-  M. Armand, F. Endres, D. R. MacFarlane, H. Ohno, and B. Scrosati. future. Nature Mater., 8:621–629, 2009.
-  K. R. J. Lovelock. Phys. Chem. Chem. Phys., 14:5071–5089, 2012.
-  Q. X. Liu, S. Zein El Abedin, and F. Endres. Surf. Coat. Technol., 201:1352–1356, 2006.
-  P. Barpanda, J.-N. Chotard, C. Delacourt, M. Reynaud, Y. Filinchuk, M. Armand, M. Deschamps, and J.-M. Tarascon. Angew. Chem., Int. Ed., 50:2526–2531, 2011.
-  V. Borgel, E. Markevich, D. Aurbach, G. Semrau, and M. Schmidt. J. Power Sources, 189(1):331–336, 2009.
-  C. Liu, Z. Yu, D. Neff, A. Zhamu, and B. Z. Jang. Nano Lett., 10(12):4863–4868, 2010.
-  C. Largeot, C. Portet, J. Chmiola, P. L. Taberna, Y. Gogotsi, and P. Simon. J. Am. Chem. Soc., 130(9):2730–2731, 2008.
-  P. Simon and Y. Gogotsi. Nature Mater., 7:845–854, 2008.
-  B. E. Conway. Kluwer, 1999.
-  A. Balducci, R. Dugas, P.-L. Taberna, P. Simon, D. Plee, M. Mastragostino, and S. Passerini. J. Power Sources, 165:922–927, 2007.
-  A. Krause and A. Balducci. Electrochem. Commun., 13:814–817, 2011.
-  R. Y. Lin, P.-L. Taberna, S. Fantini, V. Presser, C. R. Perez, F. Malbosc, N. L. Rupesinghe, K. B. K. Teo, Y. Gogotsi, and P. Simon. J. Phys. Chem. Lett., 2:2396–2401, 2011.
-  J. R. Miller and P. Simon. Science, 321:651–652, 2008.
-  S. Liu, W. Liu, Y. Liu, J.-H. Lin, X. Zhou, M. J. Janik, R. H. Colby, and Q. Zhang. Polym. Int., 59(3):321–328, 2010.
-  Y. Liu, C. Lu, S. Twigg, M. Ghaffari, J. Lin, N. Winograd, and Q. M. Zhang. Scientif. Reports, 3:973, 2013.
-  D. M. Heyes and J. H. R. Clarke. J. Chem. Soc., Faraday Trans. 2, 77:1089–1100, 1981.
-  O. Lanning and P. A. Madden. J. Phys. Chem. B, 108(30):11069–11072, 2004.
-  C. Pinilla, M. G. Del Pópolo, J. Kohanoff, and R. M. Lynden-Bell. J. Phys. Chem. B, 111(18):4877–4884, 2007.
-  A. A. Kornyshev. J. Phys. Chem. B, 111(20):5545–5557, 2007.
-  M.M. Islam, M.T. Alam, and T. Ohsaka. J. Phys. Chem. C, 112(42):16568–16574, 2008.
-  J. Kruempelmann, C. R. Mariappan, C. Shober, and B. Roling. Phys. Rev. B, 82:224203, 2010.
-  M. V. Fedorov and A. A. Kornyshev. J. Phys. Chem. B, 112(38):11868–11872, 2008.
-  M. V. Fedorov, N. Georgi, and A. A. Kornyshev. Electrochem. Commun., 12(2):296–299, 2010.
-  N. Georgi, A. A. Kornyshev, and M. V. Fedorov. J. Electroanal. Chem., 649(1–2):261–267, 2010.
-  J. Vatamanu, O. Borodin, and G.D. Smith. Phys. Chem. Chem. Phys., 12:170–182, 2010.
-  M. Pounds, S. Tazi, M. Salanne, and P. A. Madden. J. Phys.: Condens. Matter, 21(42):424109, 2009.
-  S. Tazi, M. Salanne, C. Simon, P. Turq, M. Pounds, and P. A. Madden. J. Phys. Chem. B, 114(25):8453–8459, 2010.
-  Y.-Z. Su, Y.-C. Fu, J.-W. Yan, Z.-B. Chen, and B.-W. Mao. Angew. Chem., Int. Ed., 48:5148–5151, 2009.
-  Y.-Z. Su, Y.-C. Fu, Y.-M. Wei, J.-W. Yan, and B.-W. Mao. ChemPhysChem, 11:2764–2778, 2010.
-  J. Vatamanu, L. Cao, O. Borodin, D. Bedrov, and G. D. Smith. J. Phys. Chem. Lett., 2:2267–2272, 2011.
-  J. Vatamanu, O. Borodin, D. Bedrov, and G. D. Smith. J. Phys. Chem. C, 116(14):7940–7951, 2012.
-  L. Xing, J. Vatamanu, G. D. Smith, and D. Bedrov. J. Phys. Chem. Lett., 3:1124–1129, 2012.
-  L. Yang, B. H. Fishbine, A. Migliori, and L. R. Pratt. J. Am. Chem. Soc., 131(34):12373–12376, 2009.
-  Y. Shim and H. J. Kim. ACS Nano, 4(4):2345–2355, 2010.
-  S. Kondrat and A. A. Kornyshev. J. Phys.: Condens. Matter, 23:022201, 2011.
-  S. Kondrat, N. Georgi, M. V. Fedorov, and A. A. Kornyshev. Phys. Chem. Chem. Phys., 13:11359–11366, 2011.
-  C. Merlet, B. Rotenberg, P. A. Madden, P.-L. Taberna, P. Simon, Y. Gogotsi, and M. Salanne. Nature Mater., 11:306–310, 2012.
-  P. Wu, J. S. Huang, V. Meunier, B. G. Sumpter, and R. Qiao. ACS Nano, 5(11):9044–9051, 2011.
-  G. Feng and P. T. Cummings. J. Phys. Chem. Lett., 2(22):2859–2864, 2011.
-  D.-E. Jiang, Z. Jin, and J. Wu. Nano Lett., 11(12):5373–5377, 2011.
-  R. M. Lynden-Bell, M. G. Del Pópolo, T. G. A. Youngs, J. Kohanoff, C. G. Hanke, J. B. Harper, and C. C. Pinilla. Acc. Chem. Res., 40:1138–1145, 2007.
-  J. N. Canongia Lopes, J. Deschamps, and A. A. H. Pádua. J. Phys. Chem. B, 108:2038, 2004.
-  O. Borodin. J. Phys. Chem. B, 113(33):11463–11478, 2009.
-  F. Dommert, K. Wendler, R. Berger, L. Delle Site, and C. Holm. ChemPhysChem, 13(7):1625–1637, 2012.
-  J. N. Canongia Lopes and A. A. H. Pádua. Theor. Chem. Acc., 131(3):1129, 2012.
-  Y. T. Wang, S. L. Feng, and G. A. Voth. J. Chem. Theory Comput., 5(4):1091–1098, 2009.
-  D. Roy and M. Maroncelli. J. Phys. Chem. B, 114(39):12629–12631, 2010.
-  J.I. Siepmann and M. Sprik. J. Chem. Phys., 102(1):511–524, 1995.
-  S. K. Reed, O. J. Lanning, and P. A. Madden. J. Chem. Phys., 126(8):084704, 2007.
-  T. R. Gingrich and M. Wilson. Chem. Phys. Lett., 500(1–3):178–183, 2010.
-  P. Wu, J. Huang, V. Meunier, B. G. Sumpter, and R. Qiao. J. Phys. Chem. Lett., 3(13):1732–1737, 2012.
-  A. V. Raghunathan and N. R. Aluru. Phys. Rev. E, 76:011202, 2007.
-  M. K. Petersen, R. Kumar, H. S. White, and G. A. Voth. J. Phys. Chem. C, 116:4903–4912, 2012.
-  A. S. Pensado and A. A. H. Pádua. Angew. Chem., Int. Ed., 50:8683–8687, 2011.
-  A. C. F. Mendonça, P. Malfreyt, and A. A. H. Pádua. J. Chem. Theory Comput., 8(9):3348–3355, 2012.
-  A. C. F. Mendonça, A. A. H. Pádua, and P. Malfreyt. J. Chem. Theory Comput., 9:1600–1610, 2013.
-  C. Merlet, C. Péan, B. Rotenberg, P. A. Madden, P. Simon, and M. Salanne. J. Phys. Chem. Lett., 4:264–268, 2013.
-  J. Vatamanu, O. Borodin, and G. D. Smith. J. Phys. Chem. B, 115(12):3073–3084, 2011.
-  R. Atkin, N. Borisenko, M. Drüschler, S. Zein El Abedin, F. Endres, R. Hayes, B. Huber, and B. Roling. Phys. Chem. Chem. Phys., 13:6849–6857, 2011.
-  S. Perkin, L. Crowhurst, H. Niedermeyer, T. Welton, A. M. Smith, and N. N. Gosvami. Chem. Commun., 47:6572–6574, 2011.
-  C. Merlet, M. Salanne, B. Rotenberg, and P. A. Madden. J. Phys. Chem. C, 115:16613–16618, 2011.
-  S.A. Kislenko, I.S. Samoylov, and R.H. Amirov. Phys. Chem. Chem. Phys., 11:5584–5590, 2009.
-  G. Feng, J. S. Zhang, and R. Qiao. J. Phys. Chem. C, 113(11):4549–1559, 2009.
-  C. Merlet, M. Salanne, and B. Rotenberg. J. Phys. Chem. C, 116(14):7687–7693, 2012.
-  M. V. Fedorov and A. A. Kornyshev. Electrochim. Acta, 53:6835–6840, 2008.
-  G. D. Smith, O. Borodin, S. P. Russo, R. J. Rees, and A. F. Hollenkamp. Phys. Chem. Chem. Phys., 11:9884–9897, 2009.
-  M. Mezger, H. Schröder, H. Reichert, S. Schramm, J. S. Okasinski, S. Schröder, V. Honkimäki, M. Deutsch, B. M. Ocko, J. Ralston, M. Rohwerder, M. Stratmann, and H. Dosch. Science, 322(5900):424–428, 2008.
-  R. Atkin and G.G. Warr. J. Phys. Chem. C, 111(13):5162–5168, 2007.
-  E. Hayes, G. G. Warr, and R. Atkin. Phys. Chem. Chem. Phys., 12:1709–1723, 2010.
-  S. Perkin, T. Albrecht, and J. Klein. Phys. Chem. Chem. Phys., 12(6):1243–1247, 2010.
-  S. Perkin. Phys. Chem. Chem. Phys., 14(15):5052–5062, 2012.
-  A. M. Smith, K. R. J. Lovelock, N. N. Gosvami, P. Licence, A. Dolan, T. Welton, and S. Perkin. J. Phys. Chem. Lett., 4:378–382, 2013.
-  S. Baldelli. J. Phys. Chem. B, 109(27):13049–13051, 2005.
-  S. Baldelli. Acc. Chem. Res., 41(3):421–431, 2008.
-  S. Baldelli. J. Phys. Chem. Lett., 4:244–252, 2013.
-  A. Morita and T. Ishiyama. Phys. Chem. Chem. Phys., 10:5801–5816, 2008.
-  M. Sulpizi, M. Salanne, M. Sprik, and M.-P. Gaigeot. J. Phys. Chem. Lett., 4:83–87, 2013.
-  G.-B. Pan and W. Freyland. Chem. Phys. Lett., 427(1–3):96–100, 2006.
-  G.-B. Pan and W. Freyland. Phys. Chem. Chem. Phys., 9:3286–3290, 2007.
-  M.A. Pounds. PhD thesis, University of Edinburgh, 2009.
-  Q. Dou, M. L. Sha, and G. Z. Wu. J. Phys.: Condens. Matter, 23(17):175001, 2011.
-  R. S. Payal and S. Balasubramanian. ChemPhysChem, 13(7):1764–1771, 2012.
-  K. Shimizu, A. Pensado, P. Malfreyt, A. A. H. Pádua, and J. N. Canongia Lopes. Faraday Discuss., 154:155–169, 2012.
-  X. Si, S. Li, Y. Wang, S. Ye, and T. Yan. ChemPhysChem, 13(7):1671–1676, 2012.
-  R. M. Lynden-Bell, A. I. Frolov, and M. V. Fedorov. Phys. Chem. Chem. Phys., 14:2693, 2012.
-  G. Feng, J. Huang, B. G. Sumpter, V. Meunier, and R. Qiao. Phys. Chem. Chem. Phys., 13:14723–14734, 2011.
-  M. Z. Bazant, B. D. Storey, and A. A. Kornyshev. Phys. Rev. Lett., 106:046102, 2011.
-  G. Feng, D. Jiang, and P. T. Cummings. J. Chem. Theory Comput., 8(3):1058–1063, 2012.
-  J. Vatamanu, O. Borodin, and G. D. Smith. J. Am. Chem. Soc., 132:14825–14833, 2010.
-  M. T. Alam, M. M. Islam, T. Okajima, and T. Ohsaka. Electrochem. Commun., 9(9):2370–2374, 2007.
-  M.M. Islam, M.T. Alam, T. Okajima, and T. Ohsaka. J. Phys. Chem. C, 113:3386–3389, 2009.
-  M. T. Alam, M. M. Islam, T. Okajima, and T. Ohsaka. J. Phys. Chem. C, 113(16):6596–6601, 2009.
-  F. Silva, C. Gomes, M. Figueiredo, R. Costa, A. Martins, and C. M. Pereira. J. Electroanal. Chem., 622(2):153–160, 2008.
-  V. Lockett, R. Sedev, J. Ralston, M. Horne, and T. Rodopoulos. J. Phys. Chem. C, 112(19):7486–7495, 2008.
-  M. Drüschler, B. Huber, S. Passerini, and B. Roling. J. Phys. Chem. C, 114(8):3614–3617, 2010.
-  V. Lockett, M. Horne, R. Sedev, T. Rodopoulos, and J. Ralston. Phys. Chem. Chem. Phys., 12:12499–12512, 2010.
-  M. Drüschler, B. Huber, and B. Roling. J. Phys. Chem. C, 115(14):6802–6808, 2011.
-  B. Roling, M. Drüschler, and B. Huber. Faraday Discuss., 154:303–311, 2012.
-  M. Druschler, N. Borisenko, J. Wallauer, C. Winter, B. Huber, F. Endres, and B. Roling. Phys. Chem. Chem. Phys., 14:5090–5099, 2012.
-  J. Cheng and M. Sprik. Phys. Chem. Chem. Phys., 14:11245–11267, 2012.
-  S. Zahn, J. Thar, and B. Kirchner. J. Chem. Phys., 132:124506, 2010.
-  K. Wendler, M. Brehm, F. Malberg, B. Kirchner, and L. Delle Site. J. Chem. Theory Comput., 8(5):1570–1579, 2012.
-  M. Salanne, L. J. A. Siqueira, A. P. Seitsonen, P. A. Madden, and B. Kirchner. Faraday Discuss., 154:171–188, 2012.
-  R. van Roij. In S. Dean, J. Dobnikar, A. Naji, and R. Podgornik, editors, Proceedings of the CECAM Workshop “New challenges in Electrostatics of Soft and Disordered Matter”. Pan Stanford, 2013.
-  A. Al-Zubaidi, T. Inoue, T. Matsuhita, Y. Ishii, and S. Kawasaki. Phys. Chem. Chem. Phys., 14:16055–16061, 2012.
-  L. G. Bettini, M. Galluzzi, A. Podestá, P. Milani, and P. Piseri. Carbon, 59:212–220, 2013.
-  Y. Gogotsi, A. Nikitin, H. Ye, W. Zhou, J. E. Fischer, B. Yi, H. C. Foley, and M. W. Barsoum. Nature Mater., 2:591–594, 2003.
-  G. A. Feng, R. Qiao, J. S. Huang, S. Dai, B. G. Sumpter, and V. Meunier. Phys. Chem. Chem. Phys., 13(3):1152–1161, 2011.
-  G. Feng, S. Li, J. S. Atchison, V. Presser, and P. T. Cummings. J. Phys. Chem. C, 117:9178–9186, 2013.
-  R. Dash, J. Chmiola, G. Yushin, Y. Gogotsi, G. Laudisio, J. Singer, J. E. Fisher, and S. Kucheyev. Carbon, 44(12):2489–2497, 2006.
-  J. Chmiola, G. Yushin, Y. Gogotsi, C. Portet, P. Simon, and P. L. Taberna. Science, 313:1760–1763, 2006.
-  J. Chmiola, C. Largeot, P.-L. Taberna, P. Simon, and Y. Gogotsi. Angew. Chem., Int. Ed., 47:3392–3395, 2008.
-  R. Lin, P. Huang, J. Segalini, C. Largeot, P. L. Taberna, J. Chmiola, Y. Gogotsi, and P. Simon. Electrochim. Acta, 54(27):7025–7032, 2009.
-  J. C. Palmer, A. Llobet, S.-H. Yeon, J. E. Fisher, Y. Shi, Y. Gogotsi, and K. E. Gubbins. Carbon, 48:1116–1123, 2010.
-  C. L. Bishop and M. Wilson. J. Phys.: Condens. Matter, 21(11):115301, 2009.
-  O. Pizio, S. Sokolowski, and Z. Sokolowska. J. Chem. Phys., 137(23):234705, 2012.
-  L. Xing, J. Vatamanu, O. Borodin, and D. Bedrov. J. Phys. Chem. Lett., 4:132–140, 2013.
-  S. Kondrat, C. R. Perez, V. Presser, Y. Gogotsi, and A. A. Kornyshev. Energy Environ. Sci., 5(4):6474–6479, 2012.
-  C. Merlet, M. Salanne, B. Rotenberg, and P. A. Madden. Electrochim. Acta, 101:262–271, 2013.
-  L. B. Bhuiyan and S. Lamperski. Mol. Phys., 111:807–815, 2013.
-  L. Pastewka, T. T. Järvi, L. Mayrhofer, and M. Moseler. Phys. Rev. B, 83(16):165418, 2011.
-  S. Millefiorini, A. H. Tkaczyk, R. Sedev, J. Efthimiadis, and J. Ralston. J. Am. Chem. Soc., 128:3098, 2006.
-  S. K. Reed, P. A. Madden, and A. Papadopoulos. J. Chem. Phys., 128(12):124701, 2008.