The electrochemical double layer at the graphene/aqueous electrolyte interface: what we can learn from simulations, experiments, and theory

The physical-chemistry of the graphene/aqueous–electrolyte interface underpins the operational conditions of a wide range of devices. Despite its importance, this interface is poorly understood due to the challenges faced in its experimental characterization and the difficulty of developing models that encompass its full physics. In this review we first summarize the classical theory of the electrochemical double layer, with the aim of defining a universal nomenclature to link experiments and simulations within a single unified framework. We then present the most recent experimental, theoretical and computational data and discuss how they compare with standard theory. The review ends with some remarks about how to compare simulations and experimental data and how this technology might evolve in the future.


Introduction
Nanotechnology has promised a revolution in many industrial sectors that is slowly becoming reality, but full realisation of its potential requires deeper understanding of the physics at the nanoscale. Carbon-based nanomaterials and related technologies are at the forefront of this revolution. In fact, the unique variety of thermodynamically stable configurations that carbon shows at the nanoscale (e.g. carbon-nanotubes, graphene and fullerene) allow for tailored manipulation of its mechanical, optical and electronic properties. As a consequence, carbon materials are nowadays employed in a broad range of applications including in medicine as diagnostic as well as therapeutic tools, 1 in energy storage and conversion, 2 in electronics 3 and in environmental science as sorbents, filters or sensors. 4 In the majority of these applications, the carbonaceous material is in contact with a water-based electrolyte and the functionality of the devices is dictated by physico-chemical phenomena taking place at this interface. Despite its vital importance in many areas of science, this solid/liquid interface is however not fully understood at the atomic level and in particular, the mechanism of ion adsorption/desorption onto the surface is still unclear and highly debated in the scientific community. [5][6][7] This lack of knowledge hinders the development of a wide range of technologies that rely on this physical phenomenon such as supercapacitors but also electrochemical sensors, rechargeable batteries and many optical devices. Consequently, a surprisingly large number of fundamental questions about this interface remain unanswered including: (i) what is the atomic-level structure of the electrical double layer? (ii) what are the competing interactions that drive the ions' adsorption? (iii) how does the application of an external potential affect the ions distribution and dynamics at the interface? (iv) how does confinement of an electrolyte affect the ions' dynamics and the solution electrochemical properties? The answers to these questions lead directly more technological question such as: what is the role of the ion/surface interaction on the performance of electrical double layer devices?
The emergence of graphene as a new material had created the golden opportunity to answer many of the questions posed above, since its theoretical atomic scale smoothness and single atom 2D geometrical nature make it not only the perfect substrate to perform highly controlled experiments, but also the perfect material to carry out molecular simulations to compare with experimental results. However, the lack of complexity of the experimental system is only an apparent one. Graphene often contains impurities (not only affecting its surface chemistry, but also compromising its presumed atomistic smoothness) and/or structural defects (with consequences for its electronic structure). Moreover, simulations are plagued by well-known problems associated with an uncertainty of the model (i.e. force field), limitations of fully classical simulations to model charge redistribution or charge transfer and the spatial and temporal limitations associated with the scale of quantum mechanical simulations.
The literature on graphene/electrolyte interface is vast especially concerning simulations (classical and quantum mechanical) with some valuable other review papers published fairly recently on the topic (see for example ref. 8). Here we focus solely on the work done to characterize the interface from an electrochemical point of view and the efforts that have been made to compare simulation and experimental data. We only report on water-based electrolytes with specific reference to electrolytes containing alkali metals; these have been used as model systems in experimental and theoretical studies due to the well-established fundamental theory of the thermodynamics of the electrochemical double layer.
The review starts with a summary of the classical theory of the electrochemical double layer (EDL). Despite this being textbook knowledge, we use this initial introduction to define a universal nomenclature and to link experiments and simulations within a single unified framework to make clear and unequivocal the link between the two throughout the whole review. The subsequent sections review the computational methods developed to model the interface (neutral and electrified) and the experimental electrochemical techniques for the study of the EDL, followed by a summary of the available capacitance data (experimental or calculated) and proposed ion adsorption mechanisms. We then move on to consider the effect of confinement on the EDL structure and ions' dynamics. The review ends with some remarks about how to compare simulations and experimental data and the next challenges in the field.

Classical theories of the electrochemical double layer
The term interface is used to describe the imaginary plane that separates two immiscible physical phases being in contact with each other. The thickness of such a physical interface is strictly defined to be no more than that of a complete atomic monolayer (usually on the order of a few Å) and hence the interface is treated as a two-dimensional object. If the interface extends deeper into the tangential individual phases (i.e., from at least two monolayers up to the mm scale), the term interfacial region is used instead. In both cases, the properties of the new phase formed in the vicinity of the geometrical boundary of the individual phases, are different (often significantly) compared to those of the homogeneous phases. When the individual phases are electrical conductors and at least one of them is also an ionic conductor, the interface/interfacial region is defined as an electrochemical interface/interfacial region. Interfaces of this type are formed when a metal, semiconductor, or other types of (semi)conductive materials such as graphite, conductive polymers, and various ceramics, are in contact with an ionic conductor, e.g., an electrolytic solution, solid electrolytes, ionic melts, membranes, etc.
The formation of an electrochemical interface is known to be accompanied by a spatial separation of charge, which leads to a local potential difference between the electrode surface and the electrolyte solution and perturbs the activities of the various components in the vicinity of the interface. Considering the extremely small thickness of the interface, the potential gradient arising exceeds 10 7 V cm À1 . 9 The substantial strength of this electric field is not attributed to the amount of accumulated charge at the interface (most often being on the order of a few mC cm À2 ), but to its distribution and consequently the unique structure of the electrochemical interface or, as frequently reported in the literature, the EDL.
Several phenomenological models of the EDL have been developed since the middle of the nineteenth century. The fundamental basis of these models was established from macroscopic measurements under equilibrium conditions, involving mainly the interfacial tension (electrocapillary curves; see below) and double layer capacitance. 10 The change of these parameters with the applied potential and the activities of various species in the electrolyte was monitored. The suggested theories developed to describe the system were based on a thermodynamic treatment. The general approach relies on the concept of the ideally polarizable electrode, i.e., an electrode where no charge transfer processes occur, introduced by Koenig. 11 In this case the system obeys Gibbs' equation of electrocapillarity in the form: 10,12 Here, g is the interfacial tension (the reversible work of surface formation), s the density of surface charge stored in the EDL, j the Galvani potential, G i the surface excess of the ith component in the solution and m i its chemical potential. Considering a constant chemical composition with respect to the ith component in the solution, i.e., dm i = 0, eqn (1) gives: Eqn (2) was presented by Lippmann [10][11][12] and is the basis of the thermodynamics of interfacial phenomena. Differentiating Lippmann's equation with respect to j, leads to the derivation of the formula for the specific differential capacitance of the electrode: The above equation reveals that following thermodynamic considerations we can conclude that the electrochemical interface exhibits electrical characteristics similar to those of an electric capacitor. This result is of high significance for the experimental study of the structure of the EDL, since based on the capacitance relation for a parallel plate capacitor (C = e/4pd, where e, d are the relative permittivity of the dielectric and the distance between the plates respectively), C diff can be used to provide insights into both the nature of the materials (through e) as well as the geometric characteristics of the capacitor (via d). Furthermore, eqn (3) also discloses the relationship between the electrical (charge, potential, capacitance) and energetic (interfacial tension) parameters of the interface. As mentioned previously, the dependence of g on potential has been studied in parallel to that of C towards the investigation of the structural characteristics of EDL. This was made possible by the fact that the earliest quantitative experimental studies of the electrode/electrolyte interface were primarily focused on mercury, due to the feasibility of producing a contaminant-free and reproducible surface as well as its relatively expanded, compared to other metals, potential window of stability in aqueous solvents. 12 From eqn (2) and (3) it can be inferred that for s = 0, and because C cannot be negative, the second derivative of g with respect to j is negative. This means that the change of g with j exhibits a maximum. Such plots are referred to as electrocapillary curves and their maximum value corresponds to the potential of zero charge (PZC) at the interface. The latter represents the potential at which the net surface charge at the EDL is zero and is a key parameter of the electrochemical interface as it is a property of both the electrode (through its electron work function) and the detailed composition of the solution phase. 13,14 At potentials more negative or more positive than PZC, g decreases as a result of the increased presence of cations and anions, respectively. By definition, at the PZC, C will attain its minimum value.
Based on the above thermodynamic analysis, the first and simplest model introducing the concept of the formation of a double layer at the electrochemical interface and thus describing the structure of the EDL, was formulated by Helmholtz. 9,10,12 It was assumed that no neutral species other than solvent dipoles are present and adsorption of electrolyte ions occurs solely due to electrostatic interactions, i.e., non-specific adsorption. Helmholtz envisaged that the total charge on the electrode is balanced by a monolayer of ions of opposite charge located adjacent to the surface, referred to as the Helmholtz layer. In this respect, the potential gradient is linear and steep over the thickness of the monolayer. For the region extended outside the Helmholtz layer, the model states that both its potential and composition equal those of the bulk solution, with equal concentrations of anions and cations moving randomly in the solution. In this sense, the EDL represents a nano-condenser in which the two charged planes containing the excess electrons (or holes) and the counterions are separated by a solvent layer of thickness, d H . The ratio of the electronic charge density at the electrode and the potential difference at the interface (i.e., the total potential difference required to place this charge at the interface), Dj i , gives the integral Helmholtz double layer capacitance per unit area (also known as simply integral capacitance): 10,15 In the above equation j H = Dj i À j pzc , where j pzc corresponds to the PZC as defined previously and Dj i = E + Dj ref , with E and Dj ref denoting the electrode potential and the potential drop at the reference electrode respectively. The integral capacitance is essentially an average of the differential capacitance (see right below) over the potential range from the applied potential to j pzc . Note that instead of the Galvani potential j, used in the thermodynamic analysis presented above, herein we call forth the quantity of the electrode potential which is defined as the voltage difference between the working electrode and the reference electrode. The latter is immersed in the same solution with the working electrode and has a stable, well-defined electrochemical potential. The reason for substituting Galvani potential with electrode potential stems from the fact that we aim to provide insights into the interfacial properties by comparing data for different electrolyte concentrations, in contrast to the previously made assumption of a fixed electrolyte composition. By differentiating surface charge with respect to j H , we obtain the differential Helmholtz double layer capacitance, C H : As already mentioned above, such a structure is equivalent to a parallel plate capacitor and thus C H is simplified to, where e 0 is the permittivity of free space, e H the relative permittivity of the electrolyte inside the Helmholtz layer, 16 A the surface area of the electrode and d H the distance of the closest approach of the centres of non-absorbing ionic species to the electrode surface, which is often approximated by the Debye length. In the case of specific adsorption, the characteristic length is often approximated by the size of the solvated ions instead. The Helmholtz model, described by eqn (4)-(6), has been found to qualitatively follow the experimental data for relatively concentrated solutions (usually higher than 0.1 M) in the absence of ion-specific adsorption, i.e., surface inactive electrolytes, however it exhibits two fundamental drawbacks: (i) the apparent independence of C H on the applied potential, which has been disproved experimentally for various systems and (ii) its inability to explain the dependence of C H (or the electrode charge) on the electrolyte concentration, which is observed experimentally in sufficiently dilute solutions. 17,18 In more detail, despite the fact that the Helmholtz model has been a cornerstone towards the foundation of the EDL theory, it does not consider three important physical factors, which are the applied potential, the interplay between the randomized thermal motion of the ions and that imposed by the polarity of the electrode (being more pronounced at low electrolyte concentrations) as well as the interaction between solvent molecules' dipole moments and the electrode. Gouy and Chapman postulated that even though the charge at the electrode is strictly located at the surface, this might not be the case for the solution side. In dilute electrolytes, one must consider the low density of charge carriers in the solution phase. Under such conditions, the excess charge required to counterbalance that confined to the surface of the electrode will be distributed within a larger region in the solution. This newly introduced layer will have a composition different from the bulk and it will exhibit characteristics of a diffuse layer. The potential within the diffuse layer is expected to drop smoothly from the surface of the electrode (where electrostatic forces overcome thermal processes and therefore the greatest excess charge density in the solution phase is observed) to the bulk, with a decreasing gradient across the layer. In this respect, they developed a model of the electrochemical interface adopting the following assumptions: 19,20 (i) the electrode is a perfect conductor with an evenly distributed charge on its surface, (ii) the solvent is treated as a uniform dielectric continuum characterized by e r and (iii) the ions are considered as point particles whose distribution is determined by the Poisson-Boltzmann equation. Adopting this approach for the case of a symmetrical electrolyte, the expression for the Gouy-Chapman differential capacitance per unit area is derived to be, here z is the signed charge on the electrolyte ion, e 0 the elementary charge (1.602 Â 10 À19 C), k B the Boltzmann constant (1.381 Â 10 À23 J K À1 ), T the absolute temperature (K) and j GC the potential drop within the diffuse layer. The latter can be defined in a general form as the difference between the potential at the plane of closest approach and at some point in the bulk solution. The thickness of the diffuse layer is given by the Debye length of the electrolyte, L D , which can be determined via the following equation: 10 From the above formula it follows that the larger the ionic strength of the electrolyte in the bulk solution, I 0 , the smaller the thickness of the diffuse layer. For a given value of j GC , this also results in an increase in the Gouy-Chapman capacitance (through eqn (7)) due to the compression of the diffuse layer. Furthermore, since at j pzc the surface is by definition uncharged, j GC equals the potential in the bulk solution. Therefore, it can be inferred from eqn (7) that (i) the dependence of C GC on potential will be in the form of a parabola with a minimum at j pzc and (ii) C GC at j pzc decreases with decreasing electrolyte concentration. At this point, we should emphasize that eqn (7), despite being widely quoted in the literature, cannot be used for the analysis of experimentally obtained capacitance data since the potential drop j GC is unknown. Despite addressing a major weakness (constant double layer capacitance with changes in potential) of the Helmholtz model, the Gouy-Chapman model, also exhibits weak points. Firstly, it only adequately follows the experimental data in dilute electrolyte solutions and is applicable for a narrow region near the PZC. Additionally, the model considers the ions as point charges, which leads to an anomalous rise in capacitance at high polarization (since the ions as point charges can approach the surface without any restriction due to physical barriers), while ion-ion interactions, which become particularly important at higher concentrations, are neglected. The last point is the assumption of a constant electrolyte permittivity in the whole region between the electrode and the bulk solution. This assumption is inaccurate since the relative permittivity of the electrolyte decreases when moving from the bulk solution to the diffuse layer, 16 and is in fact dramatically decreased inside the compact layer at the immediate vicinity of the electrode surface 12,16,21-23 as a consequence of the decrease in the rotational freedom of water dipoles near surfaces. 12,16,[21][22][23] The next step towards the establishment of a solid theory for the structure of the EDL, was made by Stern. 9,10,12,15 His suggestion was to combine the approaches of Helmholtz and Gouy-Chapman and thus consider the electrochemical interface as consisting of two layers. In the first compact layer, the charge on the solution side is located close to the electrode surface, i.e., the Helmholtz layer, and the rest of it is distributed in the solution within the limits defined by the compact layer and the bulk solution, i.e., the Gouy-Chapman diffuse layer. In this layer, the point charge approximation is now eliminated, and the ion-ion interactions are considered by taking the ion centres as not able to approach the electrode surface closer than a specified distance. This results in a linear potential drop adjacent to the electrode surface, following the Helmholtz model, and an approximately exponential decay in the region from the compact layer to the bulk, based on the Gouy-Chapman theory. A direct consequence of the above model (referred in the literature as the Gouy-Chapman-Stern model) is the separation of differential capacitance across the interface in two components, that of the compact layer, C H , and the one from the diffuse layer, C GC . Hence, in the same way as would be expected for two electrical capacitors connected in series, we can write the reciprocal of the total capacitance of the EDL as: 10,12,15 Following the principles of the two separate models, C H is independent from the potential, while C GC varies with the applied bias in a hyperbolic way. The total capacitance, C total , exhibits a mixed behaviour and is governed at each potential by the smallest capacitance in eqn (9). In dilute solutions and near the PZC we expect a response similar to a hyperbolic function, being characteristic of a diffuse capacitance. As the electrolyte concentration or the magnitude of the potential in the case of dilute solutions is increased, C GC becomes larger and its contribution gradually decreases, resulting in a constant C total as the overall response is progressively governed by C H .
In general, the model can explain the gross features of the total capacitance of real systems, although there are still a few points that need to be considered when the model is applied for the interpretation of experimental results. The main discrepancies arise from the apparent independence of C H on applied potential. A complete theory would require consideration of the following effects: (i) changes in the dielectric properties at the immediate vicinity of the electrode surface due to the extreme energy fields rise inside the compact layer, (ii) changes in the potential drop due to the presence of anionic and cationic excesses, (iii) adsorption (both specific and non-specific) and (iv) ion pairing phenomena. 10,15 An attempt for a refinement of the Stern model, was done by Grahame in a series of seminal works, the results of which are summarized in ref. 17. As described in an authoritative review article by the pioneers in interfacial physical electrochemistry, Damaskin and Petrii, 24 Grahame stated that the potential difference at the interface, Dj m,s (the indices m and s refer to the metal electrode and the solution respectively), can be expressed using the Galvani potentials approach as, Dj m,s = Dj ion + Dj dip + Dj met (10) where Dj ion is the potential drop related to the free electrode charge and the surface excesses of counter ions compensate this charge in the solution phase, Dj dip is the potential difference due to the orientation of solvent dipoles on the surface and Dj met the potential difference arising from the so-called ''electron tails'' that exceed the ionic frame of the metal electrode. The first term in the above equation consists of the contributions from the compact layer and the diffuse layer.
Taking as a reference point the PZC, we can write: 24 Differentiating eqn (11) with respect to charge we derive the following relation: In the above formula, the first three terms are defined as the ''effective dense layer capacitance'', 24 1/C d,eff , while C GC can be calculated e.g., for a symmetrical electrolyte, as a function of both the potential and the electrolyte concentration through, 15 where s is the surface charge density as defined previously, L B , the Bjerrum length given by e 2 /(4pe 0 e GC k B T) and s* is a constant equal to (1/8pL B I 0 ) 1/2 . The main advantage of Grahame's model is its semi-phenomenological character. 15 In more detail, the model does not directly accept the dependence of C H on the potential (and therefore the electrode charge), instead it implies that it is to be determined experimentally. On the other hand, C GC can be calculated using eqn (13) with all the parameters being known. Furthermore, the knowledge of C GC allows for the determination of C eff dependence on potential. Overall, it is considered that the Grahame model in combination with the Gouy-Chapman theory can provide a quantitative description of the EDL structure at the interface formed between various metal electrodes in solutions of several aqueous and non-aqueous electrolytes in the absence of specific ion adsorption.
Rather than the smooth exponential decay described by Gouy-Champan-Stern model, a theoretical description of the electrode-electrolyte interface at the atomistic level necessarily includes the fine structure and local arrangement of the atoms and molecules in the EDL (Fig. 1). This gives rise to a more complex electrostatic potential profile, V, which can be obtained from the charge density distribution r by solving the Poisson's equation along the surface normal, Here, r 2 is the Laplace operator. In Cartesian coordinates, taking z as the direction of the surface normal and averaging over the xy-plane, we can write the solution to the Poisson equation, These relations provide a framework for extracting the V from electronic structure and classical molecular dynamics (MD) simulations where the atoms and electronic wavefunctions are accounted for explicitly. In the case of the former, the r(r) is built from the coordinates of the nuclei and the probability density of the electronic wavefunction, whereas for the latter it is instead constructed from the atomic positions and their associated partial charges. From the electrostatic potential profile we can extract electrostatic potential difference across the interface, setting the Reproduced by permission of IOP Publishing Ltd. All rights reserved. Typical electrostatic potential along the z-axis calculated between two planar surfaces (for example graphite electrodes) separated by an electrolyte. The left electrode is at 1 V, whereas the right electrode is at À1 V, leading to a potential difference of 2 V between electrodes. When the potential reaches a plateau far from the surfaces, the corresponding value provides the potential of the bulk electrolyte V bulk , which can be seen as a reference electrode. The potential drops at each electrode DV + and DV À are also shown.
bulk electrostatic potential as the electrostatic potential of a reference electrode (Fig. 1). 25 In doing so, this lets us extract both the theoretical electrochemical double layer Integral capacitance, similar to eqn (4) and electrochemical double layer differential capacitance similar to eqn (5)

The EDL theory for carbon-electrolyte interfaces
The fundamental EDL theories, described in the previous section, were established considering a highly conductive solid phase. This is attributed to the fact that the early experimental studies of the EDL were restricted to metals with an initial strong emphasis on mercury and a gradual expansion of the research to various metals (including both single crystals and polycrystalline materials), following the developments and applications of ultra-high vacuum techniques of surface science. In the case of an electrode with a lower carrier density in the vicinity of the Fermi level compared to metals, e.g., a semiconductor or a semimetal, the above models need to be revised to include the appropriate charge contributions arising from the distribution of charge in the semiconducting solid side. The latter is known in the literature as the space-charge region and is a solid-state analogue of the diffuse layer developed within the solution phase. The potential and charge distribution within the space charge region can be described quantitatively by the Poisson-Boltzmann equation, in a similar way to the classical model for the diffuse layer of the EDL. This is reasonable since the spacecharge layer can be considered as a diffuse layer with electrons instead of ions as mobile carriers. In this respect, the differential capacitance of the space charge layer, C SC , is defined by: In the above equation the charge of the space charge layer, s GC , is determined by Gauss' law, using the integrated form of the Poisson-Boltzmann equation (from the bulk to the surface of the semiconductor) to produce the relation of the electric field. For a detailed mathematical description of the derivation of eqn (19), we refer the reader to Memming's textbook. 26 Eqn (19) can be written as, in which e s is the relative permittivity of the semiconductor, L D,SC the Debye length for the semiconductor and j SC the potential drop across the space charge layer. L D,SC is given by the following formula, where n i is the intrinsic carrier density of the semiconductor. For an intrinsic semiconductor, n i equals n 0 which denotes the carrier density in the bulk of the semiconductor (since the electron and hole densities are equal). It is noteworthy that as follows from eqn (20) the minimum of the space charge capacitance, located at j SC = 0, corresponds to a finite value for C SC (since cosh(0) = 1). This is because C SC is a differential capacitance defined by eqn (19) which is not zero at j SC = 0. 26 In a similar way to what has been described earlier, the overall capacitance of the semiconductor-electrolyte interface, C total,SC , in the absence of surface states, is given by the sum of the individual contributions of the capacitances of the space charge, the Helmholtz and the diffuse layers, following the rules for capacitors connected in series: [26][27][28] where the sum 1/C H + 1/C GC can be defined as the electrolyte capacitance 1/C electrolyte . As previously mentioned, semimetals also exhibit a lower carrier density compared to metals, which results in an additional potential drop in the solid.
In the context of carbon electrodes, Randin and Yeager 29,30 were the first to experimentally investigate the capacitance of the basal plane of stress-annealed pyrolytic graphite adopting the theory for the structure of the EDL at the semiconductorelectrolyte interface. In complete analogy with the semiconductors, the semi-metallic character of graphite gives rise to a space charge layer capacitance in the solid. The authors used eqn (22) to describe the overall capacitance of the graphite -aqueous NaF solutions in the electrolyte concentration range between 10 À5 and 0.9 M. According to the potential dependence of C SC as that described in eqn (20), they calculated the contribution of the space charge capacitance in graphite, C GR , by using the following relation, 29,30 In the above equations, C SC,GR (0) is the minimum capacitance of the space charge layer as defined for an intrinsic semiconductor (see eqn (20) and (21) for n i = n 0 ), e s,GR refers to the relative permittivity of graphite (i.e., e s,GR = 2.61 in the a-b plane and e s,GR = 3.28 along the c-axis 29 ), c is the electronic charge carrier density in graphite (i.e., ca. 5 Â 10 18 carriers per cm 3 31 ) and j s the potential at the surface of the electrode. The (DOS) in the neighbourhood of the Fermi level on graphite. 32 The same author, determined the dependence of the minimum of C SC,GR 2 on the DOS at the Fermi level, N 0 (cm À3 eV À1 ), by integration of the Poisson-Boltzmann equation in the direction across the interface using the boundaries j = j s at x = 0 and j = 0 at x -N,: 28,32 C SC,GR 2 = e 0 e s,GR e 0 N 0 (25) Based on Gerischer's approach 32 if the space charge capacitance can be deconvoluted from the experimentally determined overall capacitance on the basis of eqn (22), the DOS at each applied potential can be calculated. However, such calculations would require the exact knowledge of both C GC and C H . The former can be neglected for relative high electrolyte concentrations (often larger than 0.1 M), however the accurate determination of C H is rather challenging. Finally, in the presence of surface states, i.e., additional energy levels close to the Fermi level emerging due to the intrinsic properties of the material (such as disruptions in the crystal lattice) or as a consequence of the interaction with species at the external ambient (e.g., in the presence of adsorption), extra charge is stored in the solid phase. This accumulation of charge gives rise to an additional capacitance contribution in the solid and eqn (22) becomes, 26,27 where C SS is the capacitance of the surface states, which can be considered as a capacitor connected in parallel to C SC . C SS increases the capacitance contribution of the solid side and hence the values of the DOS. Furthermore, it can be shown that C SS also varies with the applied bias and in fact in equal directions with C SC . 26 In the classical theory of capacitors, the relationship between the charge and the voltage is expressed as, As such, the voltage, V, is proportional to the charge of the capacitor and corresponds solely to the voltage drop across the two plates of the device that are separated by an insulating (dielectric) layer. We outlined in Section 2.1 above how, in the simplest case, the EDL can be represented by such a capacitor, where the electrode and solution sides resemble the plates of the capacitor. Under this approximation the potential drop can be separated into two parts, one of which is the Galvani potential, j, that was introduced in eqn (1). For a conductive electrode, i.e., a metal, the Galvani potential is the dominant contribution. However, in the case that the electrodes are instead composed of a semiconducting or insulating material, then a second important contribution to the voltage emerges, which behaves as an additional capacitor in the system. The additional contribution is typically referred to as the quantum capacitance, C Q , or chemical capacitance. The origin of C Q can be understood by considering the electronic structure of the electrode, and in particular the band structure close to the Fermi level or chemical potential, m. In metallic electrodes the Density of States (DOS) is near infinite, so charging by applying either positive or negative bias has a minimal effect on the measured electrochemical potential difference. On the other hand, in a semiconductor, insulator or semimetal, the DOS close to the Fermi level is finite. In that situation, the effect of charging the electrode leads to the creation of electron or hole particles in the conduction and valence bands according to the energetic distribution of the DOS. Ultimately, the rate of change of the occupancy of the band edges is different from the Galvani potential, leading to the non-negligible additional capacitance in the series. Due to the nature of the electronic structure of low dimensional materials, the contribution of the quantum capacitance to the total capacitance becomes increasingly important for two-dimensional electrodes such as graphene [33][34][35][36][37] and transition metal dichalcogenides. 38,39 The total capacitance of these systems can be written as, On this basis, it should be recognized that developing an understanding of the role of the quantum capacitance in monolayer and few-layer graphene is of paramount importance in the expedition of carbon-based supercapacitors as this can be tied to the unexpectedly low capacitance of high specific surface area electrodes. 33 By employing a frequency response analysis and a simple band theory model where the electrode is modelled as a two-dimensional electron gas Ji et al. disentangled the contributions of the electrochemical double layer and quantum capacitance in 1 to 5 layer graphene electrodes. 33 The experimental analysis finds that (i) the C EDL is parabolic about the PZC, (ii) for increasing numbers of graphene layers in the electrode the rate of change of C EDL with applied voltage decreases and (iii) the minimum in C EDL generally decreases with increasing layers with the exception of the graphene monolayer. Ji et al. used a band theory picture, which does not account for any of the contributions of the electrolyte, to assess the contribution of the C Q to the total capacitance, where the total charge density is the difference between the electron and hole densities (in the electrodes), Here D(E) and f (E) are the energy dependent DOS and Fermi distribution function. The model predicts that following the graphene DOS, there is a linear increase in C Q upon +ve or Àve charging of the electrode, moreover, at the PZC there is a finite capacitance on the same order of magnitude as the experimentally measured value: B2-4 mF cm À2 . It should be noted that the residual capacitance at the PZC emerges due to a temperature-induced broadening, which is included in the band theory scheme. Unlike the experimental measurement, C Q is found to increase with the number of graphene layers and converges to a constant distribution by 5 layers. Thus, the low DOS at the Fermi level in graphene is found to suppress the capacitance at the PZC. As opposed to models of the graphene band structure, electronic structure theory simulations based on DFT can instead provide access to the graphene DOS. 34,35,37 For instance Zhan et al. employed DFT simulations in vacuum and using an implicit solvent model to describe the electrolyte, to perform an analogous study to Ji et al. 33,37 The DFT-computed C Q was also found to increase with the number of graphene layers included in the model (without converging), however, at variance with experimental measurements, the relative change in the capacitance with electrode charge is less steep for systems with fewer graphene layers.
An additional advantage of computing the C Q from the graphene DOS is the ability of this approach to capture the asymmetry between the positive and negative charging branches. Zhan et al. find that the C Q increases more rapidly on the positive charging branch, while Bhushan et al. observe the opposite trend for the differential of an integrated C Q . Our approach, based on DFTB-QMMD can also be employed to extract the C Q from the graphene DOS -factoring in changes due to the electrostatic potential of the electrolyte environment, we found only a minor asymmetry between the two branches of the differential capacitance with the negative branch increasing more rapidly. 40 This leads to the conclusion that strong asymmetries in the total capacitance arise from the C EDL component and the preferential adsorption and structuring of either cations or anions on the surface.
The direct relation of C Q to the graphene DOS, has raised additional interest towards the impact of lattice imperfections, such as topological defects, on C Q . Topological defects in graphene arise by the local rearrangement of atomic bonds that break the hexagonal symmetry of the 2D lattice, i.e., forming non-hexagonal carbon rings. 41,42 The simplest topological defects include disclinations (individual pentagons and heptagons) and dislocations (pairs of pentagon and heptagon). 41,42 Wood et al. 43 used DFT calculations to investigate the effect of various types of native graphene point defects on C Q , including Stone-Wales defects, and mono-and divacancies. Their calculations show that the introduction of topological defects in sufficient concentrations, significantly increases C Q . Further DFT calculations carried out by Pak et al., 44 probed the impact of individual point-like topological defects on C Q . Their model consisted of Stone-Wales, divacancies and diinterstitials at a defect density of ca. 0.8 at%. Assuming a marginal effect of these point-like defects on C H (which however might not be valid for higher concentrations and larger-scale defects), they demonstrated that the defective graphene layers tend to enhance C Q compared to pristine graphene. They attributed this finding to the introduction of quasi-localized states of varying degrees near the Fermi level, arising due to the disruption of graphene's extended p system. An additional noteworthy outcome of both studies is that the type of defects determines the potential range within which C Q increases; the charging behaviour of the electrode becomes asymmetric with respect to PZC. For example, Stone-Wales imperfections exhibit a much larger effect at higher voltage magnitudes and hence could be potentially more effective as positive electrodes in energy storage applications. In another very interesting work, Chen et al. 45 have studied the influence of Stone-Wales point-like defects of varying concentration on the total capacitance of the monolayer graphene/aqueous electrolyte interface, by combining electrochemical experiments with DFT calculations using the approaches described above. In line with the previous studies, they concluded that the presence of Stone-Wales defects can significantly improve the DOS of graphene near the Fermi level and therefore influence C Q . This was manifested in experimentally determined integral capacitance, which is profoundly increased for the defective graphene monolayer/4 M KCl (aq) interface when compared to the pristine graphene layer (note though that C H was assumed to be unaffected by the presence of topological defects and hence its contribution to the total capacitance of the interface was taken as constant). In addition, the rate of increase in C diff with E was found to rise significantly as a function of the defect concentration and the capacitance minimum in the C diff -E plots was increased up to 70% (1.7 mF cm À2 ), relative to the pristine value, for the graphene layer with the highest defect concentration.
A newly discovered type of topological defects in graphene that give rise to remarkable physical properties, is the twisted bilayer. 46 These defects arise from the weak interlayer interactions between bilayer graphene, which allow independent control of arbitrary azimuthal orientations between the 2D lattices. The most extraordinary example of the effect of such structural modifications on the electronic properties of graphene is the twist of a graphene bilayer from a Bernal (AB) stacking configuration to a 'magic angle' of ca. 1.11, which leads to the formation of electronic bands with very weak dispersion in momentum space, so-called 'flat' bands. 47 In a recent seminal work, Yu et al. 48 established the fundamental baseline for the understanding of the heterogeneous kinetics at the twisted bilayer graphene/aqueous electrolyte interface, by probing the dependence of the heterogeneous electron transfer rate on the twist angle using both electrochemical methods and theoretical calculations. This was achieved by investigating the effect of the twist angle on graphene DOS in the vicinity of the Fermi level. To evaluate C Q in the twisted bilayer graphene, they computed the electronic band structures over a range of twist angles. The results show that C Q is notably enhanced near the neutrality point, for twist angles ranging between 11 and 21. This is attributed to the increased population of the DOS near the Fermi level calculated in this range. Interestingly, within the same twist angle range, the fraction of the applied potential found to be confined in the EDL was increased (i.e., lower voltage drop in the twisted graphene bilayer compared to the electrolyte side), which in combination with the increase in C Q highlights the fact that the introduction of flat bands can promote a bulk metal-like behaviour in low-dimensional electrodes. Furthermore, examination of the real space distribution of the electronic states reveals that in the 'magic angle' of 1.11 the local DOS near charge neutrality is strongly localized at AA sites, whereas much lower values are calculated at AB/BA or SP stacking configuration sites. This finding was strongly observed in the significantly increased C Q values (by an average factor of ca. 4) over AA sites compared to those determined for the corresponding AB/BA sites.
Following the above discussion on C Q and its calculation using electronic structure theory, it becomes evident that an unambiguous and direct relation between theory and experiment is yet to be established. From the experimental point of view, the uncertainty associated with the accurate independent determination of C H (even at high electrolyte concentrations) renders its deconvolution (see eqn (28)) from the experimentally recorded total capacitance of the interface (see Section 2.4) at best highly approximate, if not unreliable. Commonly used assumptions considering the invariance of C H on electrode material (e.g., using C H values determined at arbitrarily chosen metal electrodes, such as Pt and Au) and/or electrolyte are far from realistic and need to be revisited. Consequently, obtaining C Q directly from experimental electrochemical data becomes a rather challenging task. On the other hand, the lack of a universally accepted computational approach for the theoretical calculation of C Q , free from approximations such as the frozen bands approximation, and from idealized (often nonrealistic) models, as well as from the inherent limitations in commonly adopted assumptions regarding the complete description of C H where its variation with applied potential and electrolyte identity is considered, hinders the correlation of theoretical findings with experimental data. In this respect, a complete understanding of the EDL at the graphene/aqueous electrolyte interface would require the development of new methods where theoreticians and experimentalists work in parallel.

2.
3 Computational methods to model carbon/electrolyte interface 2.3.1 Classical molecular dynamics. An atomic scale description of the graphene-electrolyte interface, including entropic and dynamical effects on the 100s of ns time scale and 100s of Ångstrom length scale can be obtained through atomistic models such as classical molecular dynamics (MD). A detailed description of the MD principles and algorithms can be found in many books. 49,50 In summary, MD is a computational approach, which simulates the physical evolution of a system of atoms/particles via numerical solutions to the Newtonian equations of motion. In this scheme the interactions that give rise to the forces, F i , exerted on each atom, i, are computed by means of interatomic potentials and molecular force fields (FF). During a simulation each of the individual particle trajectories is determined by integrating the secondorder partial differential equation, where m i , r i and F i are the mass, position and force associated with the ith particle at a given time step, t. The forces exerted on each atom are derived from the interatomic potential, U which itself arises from the various interactions between the particles in the system, Each of the terms that contribute to the total potential describes a specific interaction and has a specific functional form. The first term, U bond , is the bonding potential that accounts for the interactions between atoms bonded together.
The bonding potential has parametrized functions for the bond stretching, bond angles and torsional bond angles. The nonbonded interactions that arise between atoms within a molecule and between atoms in different molecules are grouped into electrostatic interactions, described by the potential U elec , or van der Waals type interactions described by the potential U LJ . The electrostatic potential is the Coulomb interaction between the (partial) charges of each of the atoms. Usually, for computational efficiency U elec is separated into short and long ranged interactions. The short-ranged interactions can be computed in real space, while the long-ranged interactions are instead calculated in reciprocal space, this is known as the Ewald summation method. All non-bonded interactions that are not considered to be electrostatic are treated as a van der Waals interaction in the classical FF. These essentially approximate the many-body electron exchange and correlation interactions that give rise to the dispersive, repulsive, and inductive forces. In many popular MD force-fields, van der Waals interactions are modelled by the Lennard-Jones (LJ) 12-6 potential, This is a pair-wise function of the atomic coordinates, with the tunable parameters e ij and s ij that govern the strength and position of the attractive potential well and repulsive part of the potential. The s ij and e ij describe the specific pair interactions between atoms of type i and j, while the minimum in the potential occurs at a distance of r ij ¼ 2 1 6 Á s ij between the two atoms. In the simulations of solid/electrolyte interfaces, one of the key interactions that governs the behaviour of ions and molecules at the interface is the attractive interaction induced by the redistribution of the electronic density within the surface and that gives rise to the double layer in the solution side mentioned in the Introductory part. The nature and strength of this interaction is determined by the surface polarizability and is a purely quantum mechanical effect. Over the years, methods to include this polarization mediated surface stabilization implicitly or explicitly into classical models have been developed. In the first case, the procedure involves the reparameterization of the Lennard-Jones non-bonded parameters s ij and e ij using the free energy of adsorption obtained from quantum chemical calculations. 51 Beyond an implicit description of the induced polarization captured in the pairwise non-bonded Lennard-Jones potential, several methods exist which account for the atomic polarization explicitly. Perhaps foremost of such methods are polarizable force fields (PFF), of which there are three leading flavours: (i) induced point dipoles, (ii) fluctuating charges and (iii) Drude oscillators. [52][53][54][55] Of these, in the modelling of graphitic surfaces the Drude oscillator model has found use in the computation of the induced polarization in the surface carbon atoms. [56][57][58][59] In the Drude oscillator approach, each atomic polarizability, a, is modelled by tethering a Drude particle to the nuclear coordinate. The Drude particle is massless and carries a fraction of the atomic partial charge q D ; in this way, the partial charge of the atom is given by q À q D . The Drude particle is bound to the atom via a harmonic potential with a force constant k D . This is a universal force constant for all atoms, which is selected such that the distance between any atom and its Drude particle is always less than any interatomic distances. In the absence of any external electric field the Drude particle oscillates around the position of the atomic coordinate. However, in the presence of a field, E, for instance the one that arises from the electrostatic potential of other atoms, the Drude particle instead oscillates around position r+d, where r is the On this basis, it is possible to define the induced atomic dipole The addition of the set of Drude particles to the atoms gives rise to a FF of the form, where U self is the harmonic potential 1 2 k D d j j 2 otherwise known as the self polarization, U bond is the standard intramolecular potential arising from the bonds, angles and dihedrals, U elec describes the Coulomb potential between atom pairs, atom-Drude particle pairs and Drude particle pairs and U LJ is the Lennard-Jones non-bonded potential. At each time step, each of the Drude particles should be first shifted to their local potential energy minima, which then allows for the selfconsistent field computation of the forces acting on the atoms in the Born-Oppenheimer approximation. The computational feasibility of optimizing the positions of each Drude particle at each time step is often questionable. Instead, a viable alternative, which approximates the potential energy surface, is to afford each of the Drude particles a fraction of the atomic mass and kinetic energy through extended Lagrangian mechanics. This allows for the treatment of the motion of the Drude particles on an equal footing to the atomic nuclei via an additional thermostat with temperature set to T * and the condition that T * { T. In this regime the Drude particles only carry a small kinetic energy that updates their positions in accordance with the thermalisation of the nuclei. Also within a fully classical setup, the electrostatic response of an electrode surface to variations in the electrolyte configuration can be approximated using the constant potential method (CPM) as introduced by Siepmann and Sprik 60 and further developed by Reed et al. 61 The CPM leverages the fact that, as we discussed at length in Section 2, for perfect conductors the internal electric field is zero and the relative permittivity is infinite, thus the electrostatic potential experienced inside the electrode should be constant. To this end, in a CPM simulation the partial charges of the electrode surface atoms are variable and determined selfconsistently at each time step to maintain an overall constant electrostatic potential throughout the electrode.
In practice, in the CPM the electrode charge density is typically described by a set of overlapping atom-centred Gaussian functions, These r j (r) have the following properties: their total integrated charge is q j , they are finite, with a width Z and are normalized by the factor A ¼ Z 3 p 3 2. The charge density arising from the electrolyte atoms is described by a set of point charges r i (r) = q i d(r À r i ), yielding a functional form for the Coulomb energy based on the system total charge density (electrode + electrolyte), Within this formalism, it is possible to express the potential that is experienced by an atomic charge at the site r j in the electrode as the derivative of the Coulomb energy with respect to the atomic partial charges, By requiring that the potential experienced by each electrode atom is a user-specified constant V 0 j (and thereby enforcing a constant electrostatic potential within the electrode) the total potential energy of the system can be found by minimizing with respect to the electrode charges q j . This minimization step has an associated computational cost reducing the scope of MD simulations with respect to the constant surface charge methods described above. Yet, the selfconsistent determination of the set of q j , in direct response to the specific configuration of the electrolyte atoms at the interface facilitates simulations on the nm and ns length and time scales that include the effects of long-ranged redistribution of metallic electrode partial charges. The importance of these effects on the structural and dynamical properties of the system at hand can be used to determine computational accuracyviability trade-offs. 62 In addition, in contrast to each of the other methods discussed here, which are all performed at constant surface charge, the CPM lets the user set the electrode potential analogous to the situation in a potentiostaticallycontrolled electrochemical experiment.
Recently, hybrid quantum mechanical-classical molecular dynamics (QM/MD) methods have also emerged as an alternative to PFFs and the CPM as a way to incorporate the fluctuations of the electrode surface charge density in classical FFs. 63 Unlike quantum mechanics-molecular mechanics (QMMM) approaches that target the properties of a quantum mechanical region embedded inside a classical electrostatic potential, these QM/MD schemes seek to describe the evolution of the classical system subject to the evolving quantum mechanical polarization of the surface.
QM/MD approaches directly couple electronic structure theory calculations of the electrode surfaces to the classical evolution of the electrolyte. As shown in Fig. 2, this is achieved by means of an iterative procedure that makes calls to QM and MD calculators. Firstly, for a given electrode-electrolyte configuration, the electrolyte coordinates are transformed to a set of point charges {q i }, these form the background electrostatic potential for the QM calculation. The QM calculation is performed considering only the electrode atoms explicitly, with the surface electron density distribution optimized according to the background electrostatic potential. In the second step, through post-processing and decomposition of the electron density, a set of atomic charges are computed that constitute the updated partial charges in the classical FF. This means that the polarization mediated interactions are captured inside the Coulomb potential of the MD steps, other non-bonded interactions such as van der Waals and London dispersion interactions remain treated in a pairwise manner via a Lennard-Jones type potential as described by eqn (34). It should be noted that this carries the implication that classical force-fields that parameterize the polarization mediated interactions, such as Williams et al., cannot be combined with the QMMD approach due to a double counting of this interaction. Evolving the entire system in time, subject to the QM partial charges, via MD establishes a new configuration for the step one.
QM/MD approaches benefit from several advantages, the in situ calculation of the electronic structure provides access to the electrode band structure, density of states and electron density distribution in the presence of electrostatic potential arising from electrolyte configuration. While there is no formal quantum mechanical operator for defining the atomic electronic occupations, the distribution of the atomic charges can be grounded in charge partitioning schemes within electronic structure theory. Therefore, this method does not require the electrode to be metallic, i.e. the QM/MD approach (unlike the CPM) can equally treat semi-metallic or semiconducting electrodes such as those made of bulk carbon materials and graphene. The procedure is however more computationally intensive than standard MD simulations due to the frequent need for QM calculations to update the surface charges. The increase in computational demand can be mitigated by adopting an analytical description of the surface polarization in the form of a so-called ''machine learning model''. To this end, our group has recently developed a procedure to train a multi-input multi-output neural network that treats the surface polarization as a discrete classification problem. The procedure is generic and introduces the electronic polarization into classical Molecular Dynamics (MD) force fields using a hybrid pairwisemachine learning force field. Assuming a good training set of data is available, the Neural Network Molecular Dynamics (NNMD) simulations can be deployed to obtain results comparable to those obtained by the QM/MD simulations but at computational complexity of a standard MD simulation.
2.3.2 Density functional theory. Despite the advances in classical models, high level quantum mechanical methods are still needed when detailed electronic structure information are sought. For example, in the QM/MD method, due to the use of Density Functional Tight Binding (DFTB) for the QM part of the workflow (the choice of the DFTB was dictated by computational cost), phenomena such as charge transfer, if they occur, cannot be properly captured during the simulation. Among the various quantum mechanical methods available, it is widely accepted that the workhorse in materials science, and in particular in the investigation of solid-liquid interfaces, [64][65][66] is Kohn-Sham Density Functional Theory (DFT). 67,68 DFT is underpinned by the fact that for any given system of interacting electrons which exist in an external potential arising from the coordinates of the atomic nuclei, the ground-state total energy has a functional form expressed in terms of the ground-state electronic density. This dependency of the energy on the electronic density, rather than the many-body electronic wavefunction, is advantageous because it drastically reduces the complexity of the problem to be solved. Formally, DFT is an exact formulation for the ground state energy, however, in practice solutions of the Kohn-Sham equations require the introduction of approximations since the total energy functional remains unknown. 67,68  The formulation of Kohn-Sham DFT relies on the two Hohenberg and Kohn theorems (1964), 69 the first of which legitimises the use of the electronic density as a variable by proving there is a one-to-one correspondence between the external potential and the electronic density. The second theorem provides a variational principle for the ground state energy. Thus, together these two theorems rewrite quantum mechanics not as a Schrödinger equation, but instead as an exact variational principle of the electronic density. The issue is that the functional is not defined. Kohn and Sham introduced a practical approach to solve the energy functional introducing a set of auxiliary noninteracting electrons, these exactly reproduce the ground state electronic density of the fully interacting system. This leads to the definition of a set of self-consistent single particle equations, known as the Kohn-Sham equations, 70 here the first term describes the kinetic energy operator with r 2 being the second derivative with respect to the coordinate system, V H is the Hartree potential, V xc is the exchange-correlation potential and V ext is a fixed external potential which describes the electron-Nuclei interactions. The set of f i (r) and e i are then the Kohn-Sham orbitals and energy levels. In Kohn-Sham DFT the total energy of the system is expressed as a functional of the electronic density, The Hartree potential which accounts for coulombic electronelectron interactions takes the form, while the quantum mechanical electron-electron interactions are grouped into the exchange correlation potential, whose functional form is not known, but may be expressed as the variation of the exchange-correlation functional, E xc [n], with respect to the density, As a result, to perform practical simulations within DFT, approximations to V xc must be introduced. In fact, many approximations to the V xc have been introduced over time; a popular library of exchange-correlation potentials adopted by many of the leading DFT software packages boasts B400 different functionals as of 2018. 71,72 Approximations to the V xc can be classified into families, and these families increase in complexity. The simplest form for V xc is the Local Density Approximation (LDA), whereby the value of the E xc at a given point in space depends on the value of the density E LDA xc = E LDA xc [n]. The description is improved by introducing a dependence of E xc on the gradient of the density, in the so-called Generalized Gradient Approximation (GGA) E GGA xc = E GGA xc [n,rn]. Meta GGAs (mGGAs) extend this trend, incorporating a dependence of E xc on further derivatives of the density and on the kinetic energy density E mGGA xc = E mGGA xc [n,rn,r 2 n,t] and are increasing in popularity.
A different approach to approximating V xc is to mix a fraction of Hartree-Fock (HF) exchange into the energy functional. This can be advantageous because HF theory treats the exchange interaction exactly, Energy functionals that include HF exchange are called hybrid functionals E Hyb . .] and can be subdivided into global and ranged hybrid functionals depending on whether a distant dependent mixing takes place or not.
Within various approximations of the V xc , solution of the Kohn-Sham equations yields the self-consistent electrondensity. Decomposition of the electron density, via charge partitioning schemes, [73][74][75][76][77] can provide information on the polarization of the electrode surfaces and charge-transfer between the ions, molecules and the surface. Furthermore, the electronic density of states (DOS) and band structure are typically readily available from the f i (r) and e i . Of particular importance for the calculation of graphene electrodes is the computation of the DOS from DFT, where i and k index the bands and Brillouin zone and in practice, the delta function is approximated with a Gaussian function. The DOS is of particular importance in the determination of the capacitance since this is directly relatable to the quantum capacitance as seen in eqn (29) and (30). While DFT based simulations provide insight not obtainable from classical simulations, they are significantly more computationally demanding than any of the other MD methods described above. This is due to the fact that formally DFT scales as the third power O(N 3 ) of the size of the system. In routine application, these are limited to systems containing several hundreds of atoms and therefore miss length and time scales appropriate for a wide range of thermodynamic and kinetic properties of the electrode interface.

Electrochemical techniques for the study of the EDL
The electrochemical study of the EDL at the interface between solid electrodes and various electrolyte media, has primarily relied on the investigation of the total capacitance of the interface with respect to the applied bias. Such studies are performed under equilibrium within a potential window where no faradaic reactions occur, that is the interface is considered as ideally polarizable. The individual capacitive contributions of the phases are deconvoluted and subsequently used for the determination of various thermodynamic and electronic parameters, based on the phenomenological models described earlier (see Section 2.1). The main techniques implemented for the experimental determination of the capacitance of the interface are electrochemical impedance spectroscopy (EIS), cyclic voltammetry (CV) and galvanostatic charge/discharge method (GCD).
EIS measures the total impedance of an electrochemical system by applying a periodic perturbation signal of small amplitude, superimposed on a constant DC potential pulse (input signal) and recording the current through the cell (output signal). Both signals are periodic functions of time and are presented in the frequency domain by means of a Fourier transform. This is in contrast to the classical electrochemical techniques where the measured parameters such as current, charge and potential are presented as functions of time.
In the simplest case of a purely capacitive interface (also known in the literature as a blocking electrode), the AC response of the system can be represented by the electrical analog of an R s,e À C total connection in series, where R s,e is an ohmic resistor related to the combined contributions of the electrolyte resistance and the electronic resistance of the cell components (i.e., electrical connections, electrode resistance, etc.) and C total a capacitor resembling the total capacitance of the interface. The total impedance of such a system in the frequency domain and in rectangular form is given by, 78 In the above equation, o is the applied frequency in rad s À1 and j is the imaginary number j ¼ ffiffiffiffiffiffi ffi À1 p À Á . It is obvious that the total reactance of the system is solely attributed to the capacitive element and thus C total can be easily determined at each applied frequency by the experimental values of the imaginary component of the total impedance, Z 0 0 total , through the relation, (42) in which f is the applied frequency in Hz. In the above relation, C total represents the experimentally determined total capacitance of the interface as defined in eqn (9), (12), (22) and (26), depending on the nature of the system under study. As frequency decreases from infinity to zero, the reactance of the system (given by the appropriately rearranged form of eqn (42)) changes from zero to infinity, since at f = 0, a DC current cannot circulate through the circuit because the modulus of the total impedance, |Z total |, goes to infinity. On the other hand, R s,e is independent of the applied frequency and hence the real part of the total impedance is always constant. In this respect, the complex plane plot of a blocking system exhibits a straight line parallel to the imaginary axis with the intercept on the real axis being equal to R s,e (Fig. 3a). The Bode magnitude plot presents at high frequencies a constant value of log(|Z total |) equal to log(R s,e ) and a straight line with a slope of À1 at low frequencies (Fig. 3b). The Bode phase plot shows changes of the phase angle, y, from zero at high frequencies (the response is purely attributed to the ohmic element) to À901 at low frequencies (the total impedance is dominated by the capacitor) (Fig. 3b). The Bode phase plots can be used to precisely determine the frequency region within which the capacitance dominates the overall AC response of the system (highlighted region in Fig. 3b). Subsequently, the total capacitance of the interface is calculated using eqn (42) and the values at each frequency are averaged to give hC total i at the applied potential.
In most practical cases, the capacitance of the system does not display an ideal behavior and hence cannot be represented by the simple model described above. The origin of this deviation arises from the presence of frequency dispersion effects. Typical phenomena leading to frequency dispersion are (i) the dispersion of time constants and (ii) kinetic dispersion effects. 79 The former is attributed to distributions of time constants along either the area of the electrode (2D) or along the axis normal to the electrode surface (3D), both resulting due to a non-uniformly active electrode area. In the second case, the dispersion arises because of the slow adsorption processes of ions and/or impurities (most often neutral molecules) on the surface of the electrode. In systems displaying frequency dispersion, the classical capacitor is replaced by the constant phase element (CPE) whose impedance is given by the formula, 79 where a and Y are the constant phase exponent (0 o a o 1) and pre-exponential factor (in O À1 s Àa cm À2 ) respectively. In general, for a 4 0.8, Y is related to the capacitance of the interface and its units can be rearranged to F s aÀ1 cm À2 . The real part in the above equation demonstrates the non-ideal behavior of the CPE, with the latter representing a leaky capacitor.
Orazem and co-workers have developed a graphical approach to calculate the total capacitance of the interface for systems exhibiting frequency dispersion effects. 80 Based on this approach, the value of the constant phase exponent, a, is determined by the slope of the log Z 00 total vs. log f plot and the effective capacitance, C eff , is then calculated at each frequency using the following equation: Several other approaches for the determination of the capacitance of the interface using EIS have been proposed in the literature. A few indicative examples include modelling of the system by means of equivalent circuit analysis in the presence 79,81 or absence 78 of frequency dispersion effects, the use of modified Nyquist plots based on the total reactance of the system (see eqn (42)) [82][83][84] and generalized phase element analysis. 85,86 CV is a non-steady state potentiodynamic electrochemical technique, which is a powerful method of studying the various physicochemical processes occurring at the electrochemical interface. It involves the application of a continuous triangular DC potential pulse versus time. During the experiment, the potential is swept linearly within a range where electrode reactions occur and subsequently the direction of the scan is reversed. In this way, the charge transfer reactions occurring at the electrode are probed during the first potential ramping, while the stability and nature of the electroactive species formed at the immediate vicinity of the electrode are also investigated upon changing the direction of the scan. In general, the response is presented as a plot of current (or current density, j g , i.e., current normalized per the nominal surface area of the electrode) versus the applied potential. The controlled parameters during a CV measurement are the potential limits, the direction of the scan, the potential scan rate and the number of potential cycles.
Complementary to the study of charge transfer processes, CV can be also used for the determination of the capacitance of the interface. The approach is based on the linear dependence of the capacitive current density, j c , on the potential scan rate, n, according to the formula, 63-65 where C diff and E are the specific differential capacitance of the electrode and the applied potential, respectively. Experimentally, this is realized by recording cycling voltammograms at several scan rates in a potential region within which no faradaic processes occur and determining C diff by the slope of the j c vs. n plot. 87,88 The value of j c for each scan rate is usually calculated by averaging the current values recorded in the applied potential region for the forward and/or reverse scan direction. 88 Furthermore, the charge density accumulated at the surface of the electrode, s, during each potential cycle can be easily found by integration of j c with respect to time, 89,90 s ¼ 1 2 where n is the potential cycle number, t 0 is the time required for the completion of half a potential cycle and E the potential of the electrode. The knowledge of s allows for the determination of the specific integral capacitance based on the formula, where E max and E min define the limits of the potential region used in the CV. As mentioned earlier (see Section 2.1), C int is essentially an average of C diff over the applied potential range. Finally, for applications in energy storage systems, the GCD method is most often used to calculate C int . The latter can be directly related to the overall performance of the device and hence is preferred over C diff . The method measures the potential of the electrode under the application of a constant current pulse for multiple consecutive cycles. C int is determined from the charging and discharging currents using the following formula, 91 where Dt is the time in between the variation of potential from E max to E min (or vice versa) under an imposed current, j c .
3. Review of the existing capacitive data obtained from experimental and computational studies

Experimental data on the graphene-electrolyte capacitance
Above we have outlined the classical and computational approaches to the electrical double layer, and the experimental methods used to interrogate the double layer, we now consider some of the experimental studies relating to the capacitance of the graphene/electrolyte interface, in particular. This interface has attracted considerable attention over the past decade or so, for a number of reasons. Graphene is the thinnest possible electrode material, thus from a ''device'' perspective, it should offer a route to maximise the gravimetric energy density (energy per unit mass) stored in a capacitor based on carbon electrodes. It should be noted that virtually all commercial supercapacitor devices are currently based on high surface area carbon materials, usually derived from activated carbon. The two-dimensional nature of graphene also raises interesting questions of a more fundamental nature, for example how does the capacitance of the material evolve with thickness? What effect, if any, does the adsorption of different ions have on the capacitance of graphene? Given that in any real device, graphene would be supported on a substrate, what effect does substrate identity or roughness have on the capacitance of the graphene-based structure? From a computational perspective, the smaller number of atoms in a ''slab'' of graphene also make it an attractive target. These factors, coupled with the general interest in graphene as a material, and in electrochemical energy storage as an approach to increase the uptake of renewable energy, have led to an avalanche of papers on the use of graphene in capacitive energy storage (the search term ''graphene and supercap-'' finds more than 21 000 papers on the Web of Science database). Clearly it is impossible to produce any kind of sensible review of this enormous body of literature, but a few salient points should be noted.
The majority of these works use graphene derived from solution-based chemical methods (such as the Hummers' method 92 ), meaning that the carbon material is more properly termed graphene oxide, or at least reduced graphene oxide, in most cases. Moreover, it is difficult to realise a practical device consisting solely of graphene sheets, therefore it is also very common for the graphene to be used as a composite material with (for example) conducting polymers such as polypyrrole 93 or metal oxides such as manganese oxide. 94 Of course, ternary structures can also be produced by simultaneously combining graphene with a metal oxide and an appropriate polymer. 95 Most of the aforementioned studies are ''top-down'', i.e. a relatively complex electrode material is formed, its (apparent) capacitance is measured, and some ex situ structural studies are performed in an attempt to justify the high gravimetric capacitance that is generally reported. How these capacitances relate to microscopic process at the electrode/electrolyte interface, given the complex physical and chemical structure of the material, is usually unclear. One recurring difficulty is the assessment of the electrochemical surface area of the porous electrode structure that is normally required for supercapacitor applications. Gas sorption measurements can give a guide, but are likely to overestimate the area exposed to the electrolyte, hence capacitance values are commonly quoted gravimetrically (i.e. per unit of electrode mass); where areal values are quoted, it is not always clear how the area has been determined. Some more discerning works have recognised this problem, and have offered routes to resolve it, further recognising that potentialdependent ingress of ions to smaller pores is likely to occur, i.e. the effective surface area is likely to be potential-dependent in many cases. 96 A much more promising approach, therefore, is to employ well-defined graphene structures, such as those formed by chemical vapour deposition methods, where a defined geometric area (and thickness) of the sample can be exposed to an electrolyte solution and the capacitance of the electrode (or electrodes in the case of a ''device'') is measured. Unfortunately, there is a relative paucity of such studies in the literature, however the key points from some of the more important studies of this type will be summarised here.
The most prominent early paper on the capacitance of monolayer graphene was the report from the laboratory of the late N. J. Tao, 97 who measured the capacitance of the graphene/aqueous electrolyte (NaF) and graphene/ionic liquid (1-butyl-3-methylimidazolium hexafluorophosphate) interfaces. The authors isolated the quantum capacitance term, which can be equated with C SC in the bulk material (see eqn (22)), by assuming that the solution components of the capacitance (C H at high ionic strengths) would be identical to those at the Pt/aq and Pt/ionic liquid interfaces. This is a questionable assumption, since it assumes that the organisation of the solvent molecules and ions are identical in both cases. Notwithstanding this assumption, the authors obtained experimental values for the quantum capacitance of graphene, which were of the order of 7-10 mF cm À2 , and could only be fitted to theory with the introduction of an effective charge carrier term due to charged impurities imparted by the substrate (oxidecovered Si).
An interesting follow-up to the work of Tao, came from the group of Ruoff, who reported the capacitance of monolayer CVD graphene where either both sides, or only one side, of the sample were exposed to an aqueous electrolyte (H 2 SO 4 ). 98 The results were interpreted within the same framework as developed by Tao et al., 97 although capacitance could be interpreted without recourse to the ''charged impurities'' invoked by Tao. Whether this is because of the difference in substrate, or because of a difference in sample preparation, is not clear at present. In fact, the article is more notable for what it omits: the two-sided configuration developed could provide extremely interesting insights on the screening of electrolyte by different thicknesses of graphene, or the effects of ion size on interactions ''across'' the graphene but the authors did not extend their study to these areas. The difficulties in obtaining consistent results for nominally the same CVD graphene configuration is illustrated by a report using CVD graphene with a polymer gel electrolyte, where an areal capacitance of 80 mF cm À2 is quoted. 99 The relatively low areal capacitance of individual graphene sheets, and the propensity of multiple graphene sheets to re-aggregate led to the development of more complex growth strategies where, for example, verticallyoriented graphene sheets were grown from conducting substrates. 100 These approaches, while interesting from the perspective of practical devices, led once again to difficulties in understanding the origin of the higher capacitance observed, since the geometries of the electrode structures are irregular, and the presence of graphene ''edges'' will substantially increase the capacitance. 101 The subsequent trajectory of the literature has been to focus on the preparation of increasingly complex structures (both physically and chemically), although a few papers have retained a focus on the capacitance of welldefined samples. For example, Chen et al. explored the effect of nitrogen-doping on the capacitance of (predominantly monolayer) CVD graphene in an aq. KCl solution. The doping imparted by the introduction of N atoms was found to increase the capacitance of the sample, particularly at potentials negative of the PZC and yielded a positive shift of the PZC. 45 More recently, however, there has been a renewed interest in more ''fundamental'' studies of the capacitance of graphene electrodes, with a focus on the effects of ion and substrate identity. For example, electrochemical impedance spectroscopy and quartz crystal microbalance methods have been confined to interrogate the double-layer formed with ionic liquids (either neat, or dissolved in acetonitrile) at monolayer graphene electrodes, prepared via CVD. 102 Specific adsorption of the imidazolium cation is the notable finding of this work, attributed to p-p interactions between the graphene and the cation. The ''substrate effect'' has also been explored in recent literature: Kwon et al. deposited CVD graphene on substrates with a range of hydrophobicity/hydrophilicity and found a trend of decreasing capacitance with increasing water-in-air contact angle. 103 This was attributed to change in the local water structure, as expressed through the local relative permittivity, due to the interaction of the interfacial water. Unfortunately, these experimental findings are completely at odds with an earlier molecular dynamics study, which had suggested that the capacitance of monolayer graphene would be largely invariant with substrate hydrophobicity. 104 The computational work did recognise that the differences in interfacial water ordering would be anticipated but predicted that these effects would be offset by changes in the ion distributions within the Helmholtz layer.
In summary, the complexities in the preparation and structure of many solution-processed graphene materials, including composite materials, makes it very hard to compare capacitance data between works. Factors such as sample thickness and porosity for these ''macroscopic'' electrodes are often not discussed, and even where they are, it is not straightforward to correct for these effects to allow for comparison between samples. Greater insight can be obtained by focussing on welldefined samples (prepared using either CVD or mechanical exfoliation) where the surface area exposed to the electrolyte and graphene thickness is accurately known.
As we have seen above, there is still some further work required to give convergence between experimental works on the capacitance of graphene and atomistic simulations. One important point, which is not always appreciated in the computational literature, is that any electrode/electrolyte interface has a finite range of stability, i.e. a so-called ''potential window''. Excursions above and below this potential window will lead to oxidation and reduction processes, respectively, which are often irreversible. These processes change the composition of the cell, which is undesirable for a device meant to be functioning solely as a capacitor. Furthermore, such processes could lead to catastrophic breakdown of the device because a conducting material can be irreversibly oxidised, resulting in material loss due to electrodissolution, [105][106][107][108][109] formation of a poorly conducting film in the case of a metal electrode, 110,111 or the evolution of gaseous products, e.g. CO 2 , in the case of a carbon electrode. 112,113 The potential window of each interface is governed by the electrochemical stability of the least stable component, from the electrodes, solvent and electrolytic salt. Water has a particularly limited range of stability, both oxidative and reductive: the thermodynamic limit of water's stability gives a potential window of 1.23 V, 9,114 although this is extended by kinetic factors due to the slow electrolytic breakdown of water on carbon electrodes in particular. 105,112,115 For this reason, organic-based electrolytes are used in commercial supercapacitors, with acetonitrile combined with tetraalkylammonium salts of anions such as tetrafluoroborate, being a common choice. This solvent/salt combination, coupled with activated carbon electrodes, gives an operational voltage of 2.7 V for commercial supercapacitors, 116 although research is based on improving the stability of the components to increase this voltage and thus enhance the energy density of the devices. 117 The voltages quoted above should be borne in mind for computational studies. Occasionally, substantially higher voltages are mentioned in the literature, 118 but these will lead to cell breakdown, at least for currently known electrode/electrolyte combinations and are therefore simply ''unphysical''. One interesting, but unresolved, area of research is how the electrochemical stability of graphene compares with that of the activated carbons used in current commercial devices. In principle, this factor should be easier to determine than it is with activated carbons, due to the heterogeneity of the latter.

Computational studies on the double layer structure and ions adsorption mechanisms
The graphene/electrolyte interface is also more amenable to accurate simulation, than the physically and chemically complex activated carbon materials. Simulations provide the bulk of the understanding we have currently about the atomic scale properties of the graphene/electrolyte interface. The graphene-water interface has been intensively investigated 119 in the context of both MD 56,120-126 and ab initio simulations. [127][128][129][130][131][132][133][134][135][136][137][138][139] Here we focus on the electrochemistry of the interface, which necessarily involves the physisorption and chemisorption of ions from the aqueous solutions onto the electrode surface. The addition of ions, and the inclusion of surface charge, changes and increases the complexity of the problem at hand; in MD simulations the additional ion-water and ion-surface interactions are a delicate balance, whereas in DFT simulations, the trade-off between chemical accuracy and computational viability becomes ever more important.
Chialvo and Cummings employed constant-pressure constanttemperature (NPT) MDs to characterize the graphene-electrolyte interface for cations of different ionic charge and radius and subject to confinement. 140 In their work they considered finite, free standing graphene sheets immersed in both pure SPC/E water (simple point charge model with rigid bonds) and aqueous 2 M LiCl, BaCl 2 and YCl 3 solutions. The C atoms within graphene were treated as non-polarizable Lennard-Jones spheres using Steele's parameters, 141 with C-O and C-ion interaction potentials determined via Lorentz-Berthelot combining rules.
Plots of the atom-resolved density distribution along the direction of the surface normal reveals that for each LiCl there is preferential adsorption of the cation closer to the surface than the Cl À anion. This gives rise to the characteristic double layer structure. In addition, an important feature that arises from Chialvo and Cummings' simulations is the retention of the ion hydration shell in the adsorbed configuration.
Their density plots show that oxygen atoms belonging to the SPC/E water molecules always adsorb closer to the graphene surface than the Li cations. The adsorption of the various cations in their fully hydrated configurations leads to so-called ''specific ion effects'', where the smaller Li + ion has a strong adsorption peak 0.4 nm from the surface, while larger Ba 2+ and Y 3+ ions behave differently. The intermediate Ba 2+ ion has a similar, but significantly less intense, adsorption peak at 0.55 nm from the surface and a strong adsorption at 0.8 nm. Finally, Y 3+ being the largest ion investigated has a first adsorption peak at 1.5 nm, thereby inverting the charge density distribution of the double layer by residing further away from the graphene than the Cl À counterion.
Recognizing the importance of surface polarization effects, and in particular the effect of cation-p orbital stabilization on the renormalization of the graphene-ion interaction potential Williams et al. employed DFT simulations in order to parameterize a set of Lennard-Jones parameters suitable for Li + , Na + , K + , Mg 2+ and Ca 2+ chloride salts adsorbed on a graphene layer. 51 To model the interactions between ions and graphene Williams et al. used a hexagonal graphene flake, terminated with hydrogen atoms at the edges and a single cation adsorbed directly in the centre of the flake. To account for the role of the water molecules, they employed an implicit solvent model based on the conductor-like polarizable continuum model (CPCM), with modified cavity radii optimized to match experimental hydration free energies. Inspection of the DFT adsorption free energies suggests that in the aqueous phase, all cations considered, and the chloride anion, prefer to adsorb on the surface rather than remaining in the bulk phase. In the gas phase, Williams observed a clear trend with the ion adsorption energy becoming increasingly negative with smaller ionic radii, which suggests more favourable adsorption. In the aqueous phase adsorption energies are significantly damped, with no identifiable trends in ionic radius. The divalent ions Mg 2+ and Ca 2+ have stronger adsorption energies than Li + , Na + and K + , and Cl À has a significantly less negative adsorption energy than each of the cations.
Based on the DFT adsorption profiles, Williams et al. derived a set of Lennard-Jones parameters for MD simulations by matching adsorption free energies computed from MD potential of mean force 142 calculations to the optimized DFT adsorption energy. As such, the DFT adsorption energy for each specific ion was rigorously computed by employing an implicit solvent model 143 which was modified to match the experimental free energy of hydration for each ion. 51 Subsequent MD simulations of an electroneutral infinite graphene sheets in contact with 1 M electrolytes for each of the cation chloride salts permitted inspection of the interfacial structuring of the ions. Williams et al.'s force field predicts that monovalent cations always adsorb closer to the surface than the Cl À ions and that in the adsorbed configuration, Na + and K + are partially dehydrated suggesting direct adsorption of the ion on the surface. In contrast to the monovalent ions, the divalent Mg 2+ and Ca 2+ are found to reside further from the surface than the Cl À , which contributes to them retaining the full hydration shell.
Several groups have adopted this parameter set to investigate different facets of the aqueous ion-graphene interface. [144][145][146] For instance, Dočkal et al. have examined the water and ion structural and diffusive properties within the interfacial layers. They found strong, densely packed layered structure up to 1 nm from the surface, and that the cations act as structure makers. These structure makers are found to slow the diffusion coefficients of water molecules with increasing salt concentration. 145 Within the interfacial layers where the ion concentration is highest, water molecules are most greatly slowed. Very recently Dočkal et al. also considered the cation behaviour on charged graphene sheets, 146 in summary, while the long-ranged structuring of the water layers are largely unaffected by the charged state of the electrode, locally the hydration shells of the Na + and K + ions are found to be disrupted. Through the calculation of the screening factor, counterintuitively it was found that the concentration of Cl À ions at the negative electrode increases due to the increase in concentration of positive ions (and vice versa) creating a surplus of ions. This surplus of ions can lead to a situation where water dipoles are oriented with the dipole pointing (in the conceptually opposite direction) at the surface. Finally, Finney et al. employed the NaCl parameters to investigate both surface charging and ion concentration effects on the computed electrode capacitance. 144 In particular, they also found that the Cl À ion, has an increased concertation close to the negative electrode, forming a charge compensation layer. Finney et al. were able to uncover a switch in this EDL charge screening behaviour that occurs as the concentration exceeds 0.6 M. Above 0.6 M alternating layers of positive and negative ions form in decreasing concentration until the surface charge is completely screened. Finally, they subtracted the computed double layer differential capacitance, obtained through eqn (18), from the measured specific electrode differential capacitance to gain insight on the role of the quantum capacitance. Finney et al. conclude that for graphitic electrodes close to the potential of zero charge, the relative magnitudes of the C EDL and C Q are the same; as such they point out that even at moderate concentrations the capacitance does not report directly on the electrode density of states.
Prior to their work on the behaviour of ions at the grapheneelectrolyte interface, Misra and Blankschtein developed a polarizable force field capable of capturing the self-consistent polarization of the surface C atoms based on the Drude oscillator approach. 56 In the design of their force field, Misra and Blankschtein began by considering the adsorption of monomeric water molecules onto polyaromatic hydrocarbons of increasing sizes (C 6 H 6 , C 54 H 18 & C 96 H 24 ) via ab initio simulation methods. By benchmarking against Coupled Cluster Singles Doubles and perturbative Triples (CCSD(T)) and quantum Monte Carlo approaches, they determine that the water molecule binding energy obtained at the symmetry adapted perturbation theory (SAPT0/jun-cc-pVDZ) level could be reliably extrapolated to the case of infinitely periodic graphene. Decomposing the SAPT0 total energy into its constituent contributions, Misra and Blankschtein parameterized a classical water-graphene LJ potential according to the quantum mechanical charge penetration (charge transfer), electron exchange and dispersion interactions. It is worthwhile to highlight that they removed the contributions due to the chargemultipole interactions in polyaromatic hydrocarbons, that are known to lead to under-and over-binding in periodic graphene. 63 They surveyed several different adsorption orientations for the H 2 O molecule, in spite of the fact that the SAPT0 predicted a larger binding energy for water molecules with 1 or 2 hydrogen atoms pointing towards the surface, ensuing MD simulations (both with nonpolarizable and polarizable carbon atoms) instead indicated that the stabilization brought about by maximizing the interfacial hydrogen bonding network leads to water molecules that lie parallel to the surface. This finding is consistent with the latest experimental surface sensitivespectroscopic measurements. 147 Consequently, the LJ potential was parameterized according to the adsorption profile for the parallel (tangential in their work) configuration.
The contribution to the water molecule binding energy due to polarization effects can be separated into two components: (i) the polarization of the carbon atoms in the graphene layer in response to the water molecule electric field and (ii) the induced polarization of the water molecules, Misra and Blankschtein employ the Drude oscillator model. 56 They identify two parameters, a C ; the static dipole polarizability of the C atoms and a C ; the damping parameter for electrostatic interactions between the Drude particles and the nuclei introduced by Thole. These parameters can be used to build the dipole field tensor, T, associated with two interacting dipoles, m p and m q , and therefore can be tuned to reproduce the graphene polarizability including the anisotropy between the in-and out-ofplane response of a single sheet. Here a and b run over each of the Cartesian axes, R a = a p À a q is the distance between the dipoles along the Cartesian axis a and d is the Kronecker delta. Misra and Blankschtein use the exponential dipole field tensor (which reduces to eqn (49) in the limit a -N), Here, the and a i is the dipole polarizability of the ith atom. The parameter a is the Thole damping parameter. By computing the induced polarization in a graphene cell, via the induced dipole moments of the C atoms and leveraging reported computed Berry's phase values for the out-of-plane polarizability a > = 0.837 Å 3 of a single graphene sheet and anisotropy ratio a k a ? ¼ 3:5, Misra and Blankschtein derive a C = 1.139 Å 3 and a C = 1.507 parameters for the Drude oscillators that exactly reproduce the polarizability tensor for a periodic graphene sheet. Comparison of the developed PFF with SAPT0 results revealed an underestimation in the computed induction energy of water molecules on the surface, but not with a Na + or Cl À ion. The underestimation arises in the case of water due to quantum mechanical charge penetration (transfer) effects at close distances that cannot be captured classically. Instead, a refitting of the parameters a C = 1.80 Å 3 and a C = 0.50 according to the SAPT0 induction energy yields a polarizable force field in excellent agreement with the ab initio results.
Turning to the polarization of the graphene by ions, Misra and Blankschtein applied their approach to the behaviour of SCN À ions at the graphene/electrolyte interface (Fig. 4). 58 SAPT0 calculations carried out in vacuum indicate a strong adsorption of the ion due to local polarization of the carbon atoms. Parameterization of the (Drude oscillator model) and ensuing simulations of the SCN À electrolyte/graphene interface reveals that the strong polarization of the graphene by the SCN À is significantly dampened in the presence of water. The fields generated by the ions and the water molecules are found to destructively interfere, which leads to a water-mediated screening of the ion. Interestingly, analogous calculations using a static polarization model, with a specifically parameterized LJ potential, and surface charge distributions that do not change over time, fail to capture this water screening leading to an overbinding of the ion, a clear limitation of this static polarization. This ultimately leads to polarization results being equivalent to no-polarizations results, however, in the event that the screening by the solvent is reduced (i.e. in nanochannels, ion intercalation, nanotubes, etc.) these results will diverge and a dynamical polarizability will likely be required.
Most recently, Misra and Blankstein have extended their optimized Drude Oscillator approach to survey the properties of several different anions: 59 SCN À , NO 3 À , Cl À , F À and SO 4 2À that make up the Hofmeister series. 148,149 These are ions that increase in hydration free energy moving from SCN À to SO 4 2À , which leads to a decrease in protein solubility in solutions of those ions. As with their previous work, the focus is on the ion and water-dipole induced polarization effects on the carbon atoms in the graphene surface, and the ensuing feedback which governs the ion adsorption behaviour.
In the first instance, they consider the extremes of the Hofmeister series, the highly kosmotropic anion SO 4 2À and the highly chaotropic anion SCN À . Static QM simulations, based on the SAPT0 approach, of a single ion adsorbed on the graphene surface in vacuum (note that graphene is modelled as a polycyclic aromatic hydrocarbon), indicate that both ions adsorb with a favourable negative change in the Helmholtz energy. Decomposition of the SAPT0 total energy suggests that the induced-polarization energy dominates the interaction between the ion and graphene.
In addition, in line with the expected trend that the polarization energy varies with the square of the ion charge, the SO 4 2À ion is found to adsorb with approximately 4Â the energy of the SCN À . Misra  in vacuum, and used this as a platform to investigate the effects of all of the ions in solution interacting with a periodic graphene layer. We note that these simulations are still carried out using a single ion, with the center of mass kept at fixed height above the graphene, but allowed to evolve in the in-plane directions.
In solution only the chaotropic ions SCN À and NO 3 À are found to favourably adsorb on the surface, at approximately 0.35 nm; Cl À , F À and SO 4 2À are instead all repelled from the graphene layer. Examination of the electric fields generated by the ions and water molecules moving in the vicinity of the surface suggests the formation of so-called molecular waves. These are formed from the ions and water molecules as they move over the graphene surface. The molecular waves interfere destructively to significantly reduce the polarization mediated stabilization (observed in vacuum). Consequently, Misra and Blankstein find that in solution the adsorption behaviour of the ions is instead dominated by the ion-water and water-water interactions, which explain the trends observed for adsorption of the Hofmeister series of ions. The interaction of the water molecules with the graphite and graphene dominates for small ionic concentrations, driving the wetting of the surface and thus potentially affecting the ions adsorption. Interestingly the investigations of the wetting phenomenon highlighted one of the major differences between graphene (i.e. one or few layers of one atom thick material) and graphite interfacial properties. While graphite, as any bulk material, has a well-defined water contact angle (WCA) in air, the graphene WCA changes depending on the nature of the supporting substrate and the number of graphene layers. 126 Assuming there are no airborne contaminants affecting the surface wettability, 150 experiments and simulations have indicated that the wetting of graphene depends on the wettability of the supporting substrate as the contact angle is affected by both the liquid-graphene and, due to the graphene angstrom-scale thickness, the liquid-substrate interactions. [151][152][153] In some cases, this solvent-substrate correlation across the graphene layer not only affects the WCA 154 but also the ions' adsorption. Pykal et al. employed a polarizable Drude oscillator model to inspect different graphene/electrolyte interfaces. 155 In particular, they compared differently structured interfaces where graphene was in contact with (i) electrolyte on both sides (EGE), (ii) electrolyte and pure water (EGW) on different sides, and (iii) electrolyte and vacuum (air) (EGA) on different sides (see Fig. 5).
The Drude Oscillator model they employ is adapted from the work of Ho and Striolo, 120 with isotropically polarizable graphene. Furthermore, the Drude Oscillators are fine-tuned by means of the Thole parameter. In their work, Pykal et al. avoided salt precipitation by tuning the Thole parameter of the SWM4-NDP water 156,157 to match the experimental osmotic pressures of the KF and KI solutions.
For the electrolyte, Pykal et al. investigated the behaviour of KF and KI solutions at the three differently structured interfaces (EGE, EGW and EGA). The concentration of these solutions was selected to be approximately 1.2 M; simulations carried out at 0.6 M were used to verify that there are no effects of the concentration on the ion properties at the interface. The KF and KI were chosen since the dissolved F À and I À ions exhibit qualitatively different behaviour at other interfaces As shown in (a), the electric field exerted by the SCN À ion in vacuum polarizes the carbon atoms in graphene (shown by the black solid arrow), which in turn exert back electric fields onto the SCN À ion (shown by the green dashed arrow). The induced dipole moment of carbon atom i in graphene, l ind i , is a function of the electric field exerted by the SCN À ion and the electric fields from all of the other induced dipoles in graphene. On the other hand, at the graphene/water interface shown in (b), both the SCN À ion and water molecules exert electric fields, which polarize the carbon atoms in graphene. As an example, we show the electric fields exerted by the SCN À ion and a highlighted water molecule, j, on carbon atom i in graphene (shown by black solid arrows). The polarized carbon atom i, in turn, exerts back electric fields onto the SCN À ion and the water molecule (shown by the purple dashed arrows). Because l ind i depends on both the electric field exerted by the SCN À ion and the electric fields exerted by all of the water molecules in the system, the graphene-ion and graphene-water interactions are therefore coupled. Color code: sulfur in SCN À : orange; carbon in SCN À : violet; nitrogen in SCN À : black; carbon in graphene: blue; oxygen in water: red; and hydrogen in water: white. The directions of the electric field vectors are for representation purposes; we adopt a sign convention which assumes the electric field vectors to emanate from a positive charge distribution.
(such as the electrolyte/air interface) due to the interplay between the adsorption and hydration free energies.
Across the different interfaces, the simulations of Pykal et al. suggest similar adsorption behaviour for H 2 O, K + and the two different anions. As depicted in the bulk-normalized density plots of Fig. 5, the water adsorbs with characteristic two-layer structure, while I À always resides preferentially in the interfacial region and F À in the bulk. The different behaviour of the F À and I À can be understood from the fact that the hydration free energy of the F À ion is significantly more negative than the I À and that taking a F À from the bulk phase to the interface, thence removing water molecules from its solvation shells, is a net gain in energy. 58 Both anions are found to have an adsorption peak closer to the graphene than cation K + , a result which is at odds with many of the other investigations carried out for neutral graphene-electrolyte interfaces. 51,63,158 This could stem directly from the different (and modified) water model 155,156 and specific ion parametrization 157 or more generally from the use of Drude Oscillators to model the polarization.
The fine structure of the first few layers, closest to the graphene, is found to vary between models. On the one hand, the EGE and EGW models yield similar results, with a denser packing of the water molecules at the surface. On the other hand, the EGA model, has less dense packing on top of the surface. This difference is attributed to attractive dispersion interactions between the water molecules through the graphene sheet in the case of the EGE and EGW models. Interestingly, Pykal et al. highlight that the hydrophobicity of the graphene contributes to an effective shielding of ions from the surface. Ultimately this leads them to conclude that ion-ion interactions through the graphene layer do not take place, and an independent ordering of the ions is observed for all the systems.
We developed our QMMD approach and used it to investigate the behaviour of 1 M NaCl solutions in contact with (i) hexagonal graphene flakes whose edges are hydrogenated, and (ii) free-standing infinite graphene sheets. First, we quantified the magnitude of the polarization of the graphene flake in the presence of an ion by computing the surface Mulliken charges in response to a point charge fixed at different heights above the graphene layer. We carried out this analysis in vacuum at the DFT level (6-31G/B3LYP) [159][160][161] and at the DFTB level (mio-1-1), 162 leveraging the symmetry of the graphene flake to measure the redistribution of charge amongst the rings of atoms as depicted in Fig. 6.
From this we observed a long-ranged redistribution of the Mulliken charges, which becomes stronger with the proximity of the charge to the surface. We were able to identify radial trends in the polarization state of the surface: (i) an accumulation region in the carbon rings closest to the adsorbing atom where there is an excess of charge density, (ii) a buffer region which is invariant to the presence of the ion and (iii) a depletion region defined by the rings furthest away from the adsorption site (and closest to the flake edges), which registers a loss of charge density. It is worth noting that the charge accumulation region can be as far reaching as three carbon bonds, which demonstrates the non-locality of the overall polarization gradient.
For dynamical investigation of the graphene-ion systems it is important that there is a smooth handshake between the QM and Classical parts of the simulations, in particular that the electrostatic forces acting on point charges due to the Mulliken charges are conserved when treating the electrostatic forces between the ions and the surface partial charges in the classical system. We verified from the derivative of the DFTB total energy, where Dd is the distance between a point charge and a graphene surface, and through explicit computation of the vectorial Coulomb force where Q is the charge on an ion, the set {q a } are the partial charges on the surface in the classical force field and r Qq is the vector separating Q and q a , that the force acting on a point charge/ion is conserved to within 0.02 kJ mol À1 even at unphysically close adsorption distances. In particular, at 1.0 nm distance the force obtained via both methods is À0.35 kJ mol À1 nm À1 , close to the optimal adsorption height (0.3 nm) the force from DFTB (classical Coulomb) was À27.55 (À27.54) kJ mol À1 nm À1 , while at very close distance, 0.2 nm, the DFTB (Classical) force was 87.45 (87.43) kJ mol À1 nm À1 . The method was applied to the case of 0.5 M and 1 M NaCl solutions in contact with the graphene flake, with a view to understanding the structuring of the cation and anion at the interface. Na + cations were found to reside within the water layers closest to the surface, while Cl À tends to be in the bulk water. Decomposition of the non-bonding interactions in the classical force-field shows that for Na + the attractive Coulomb interaction overlaps with the attractive part of the Lennard-Jones interaction profile, whereas for Cl À the two components always cancel out. These results suggest that cations may adsorb more favourably on graphene surfaces as has already been discussed in this section of the review. Concerning the local structure of the Na + ion at the interface, the QMMD simulations revealed that the Na + ions preferentially adsorb indirectly on the surface, that is with the full first solvation shell of surrounding water molecules. It is thought that this arises from the relatively strong hydration free energy of the Na + ion, the non-bonded interactions between the ion and the surface are not strong enough to overcome the interactions between the water molecule and the ion and thus it remains hydrated.
Finally, we considered the effect of charging the graphene surface to see how this affects the ion adsorption properties; to this end, we considered an infinite graphene sheet. Negative charging of the surface was found to increase the density of Na + ions in the water layers close to the interface, but not to induce dehydration of the ion. On the positively charged surface, the density of Cl À ions in the EDL is also found to increase, with the position of the main peak in the density profile relative to the graphene found to occur in the same place.
We further extended this work to consider different water models and a host of different ions pairs, 40 namely the TIP4P-2005 water model 163 that is known to best describe the interactions between graphene and water 126 and the Madrid-2019 ions Li + , Na + , K + , Mg 2+ and Ca 2+ . 164 We considered the local adsorption properties for each of these ions on the infinite graphene sheet under neutral and biased (charged) conditions, the density profiles are reported in Fig. 7.
Following on from our earlier work, we found similar properties for each of the ions when the surface is electroneutral, i.e. that cations preferentially adsorb closer to the surface than Cl À anions and that they retain their solvation shells adsorbing indirectly. However, owing to its comparatively low hydration free energy, upon negative charging of the graphene, the K + cation is observed to lose one water molecule and directly adsorb to the surface approximately 0.5 nm closer than the neutral case. Investigation of the dynamical properties of the various ions revealed that the direct adsorption of K + leads to a decrease in the in-plane diffusivity when the surface is charged, 40 while for the other ions the diffusivity is relatively unaffected by the surface and the charged state.
Intuitively, it is reasonable to expect that the fundamentally different adsorption mechanism and dynamical behaviour of the K + ion could lead to differences, sometimes referred to as specific ion effects, in the measured capacitance. Yasuda et al. 165 have recently identified experimentally such phenomena on monolayer graphene on atomically flat Au(111) substrates. By using electrochemical techniques, in situ Raman spectroscopy and crystal truncation rod surface diffraction measurements, they show that K + adsorption on graphene triggers electron doping with lattice expansion. The latter leads to a significant increase of the interfacial capacitance within the potential window where the adsorption process occurs. Our simulations however, revealed that the EDL contribution to the electrode integral capacitance is constant across the monovalent ions Li + , Na + and K + , a feature we attribute to the fact that the potential drop, obtained from eqn (16), involves integration over the entire double layer and so long as the concentration of ions in the EDL is constant, then the Dj could reasonably be expected to be constant also.
Colherinhas et al. made use of DFT simulations in order to understand how the adsorption of cations on the surface impacts on the electronic structure of graphene. 166 They modelled graphene, in open boundary conditions which necessitates H-termination of the graphene edges, with circumcoronene and used the BLYP GGA exchange correlation potential 167 and LANL2DZ Gaussian function basis set. Within this setup the socalled graphene exhibits molecular properties such as a finite energy gap between the Highest Occupied (HOMO) and Lowest Unoccupied (LUMO) Molecular Orbitals. Colherinhas et al. looked at the adsorption behaviours of single Li + , Na + and Mg 2+ cations in the absence of any solvent effects. They track the relative energetic alignment of the empty metal orbitals with respect to those of the graphene as a function of the ion adsorption height. An interesting outcome of their work is that at particularly short distances (o0.3 nm) the LUMO associated with the Li + and Na + ions drops below the LUMO of the circumcoronene, thereby in effect tuning the size of the electronic gap. It is posited that this effect could therefore be brought about by applying external pressure on the system. Instead, Mg 2+ does not replicate this behaviour due to a much stronger interaction with the graphene than the monovalent ions. In this case, hybridisation of the Mg 2+ and graphene molecular orbitals leads to charge transfer interactions and a 4 times stronger adsorption energy.
A similar investigation was carried out by Sangavi and co workers, who investigated the adsorption of Li + , Na + and K + cations on graphene and defective graphene structures. 168 Again, due to complications associated with charged systems, Sangavi et al. modelled their graphene in open boundary conditions, with terminating H atoms and in the absence of solvent effects. They considered a 6 Â 6 rhombohedral arrangement of the C atoms and constructed models of di-and tetravacancies as well as the well-known Stone-Wales (SW) defect. Sangavi et al. compared the results of the M06-2X hybrid metafunctional 169 with the PBE GGA, 170 using the 6-311G(d,p) and 6-31G(d) basis sets respectively. 171,172 The introduction of defects in the structure leads to the formation of Ångstrom scale pores in the hexagonal lattice, which increase in size from the SW to divalent and finally tetravalent defects; in fact defectinduced strain leads to the rippling of the tetravalent graphene sheet. It is not clear whether such ripples would be present in larger or infinite graphene models. The simulations of Sangavi et al. suggest that while the optimised adsorption distances do not vary across the graphene models, the monovalent ions, and in particular Li + , adsorb more favourably onto the di-and tetravalent defects. Li + is adsorbed most favourably due to a strong electrostatic attraction between the Li ion and defect site brought about by the purely ionic bond formed.
On the other hand, similar investigations of anions adsorbed on graphene have been carried out by Xiaozhen et al.; 173 they carried out DFT simulations on several different classes of anion spanning the Hofmeiester series. Structures were optimised with the B3LYP hybrid exchange correlation functional, 174 while the adsorption energies for each of the ions on top of graphene were computed using the oB79X-D range separated hybrid functional, 175 empirically corrected for dispersion interactions and 6-31G+(d,p) Gaussian basis set. 171,172 To model the graphene, Xiaozhen et al. considered a rectangular polyaromatic hydrocarbon with H-termination at the boundaries. A key takeaway from the work of Xiaozhen et al. is that anions with a higher ionic charge adsorb more strongly to the graphene surface due to increased mixing of the graphene and ion orbitals. Specifically, of the ions tested HPO 4 2À and SO 4 2À were found to have the most negative adsorption energy. Xiaozhen et al. have attempted to account for the role of ion hydration on the properties of the adsorbed ionic species by including one hydrating water molecule in their simulations. While, the effect of the water was to lower the adsorption energies, the trend in adsorbing ions with a higher charge more strongly was preserved. This result is in stark contrast to the classical molecular dynamics models outlined above, where higher charged ions do not adsorb strongly on the surface due to high energy barriers preventing their dehydration and exposure to the graphene surface.
In two separate works Olsson et al. 176,177 instead considered the adsorption and mobility of (atomic) Li, Na, and K, and Mg, Ca and Zn on graphene in periodic boundary conditions and using a plane waves basis set. Olsson et al. treated their graphene models with the PBE 170 exchange correlation potential, including Grimme's D3 empirical dispersion correction scheme. 178 Simulating the adsorption of atomic metals on the graphene sheet overcomes complications associated with charged simulation cells, but describes a chemical situation somewhat different to those described above. For example, in the case of the monovalent metals, 176 Olsson et al. report that Li and K exhibit preferential adsorption on the graphene surface compared with Na. This is linked to the ability of the metal to transfer charge to the graphene: For Li, Na and K the energy of the valence s-orbital relative to the system Fermi levels are 1.11, 0.27 and 0.57 eV respectively. Consequently, the Na ion that forms can most easily revert to its atomic state and is least bound to the surface. In the case of the divalent metals, Olsson et al. observe that Ca adsorption is favoured in comparison to Mg and Zn metals. 177 The adsorption of Ca is so favoured that its migration across the surface, measured by means of Nudged Elastic Band calculations, is hindered.
Leveraging the insight from electronic structure theory simulations, Zhan et al. have also explored the role of specific ion effects in the structuring and properties of the graphene/ electrolyte interface. 179 As discussed in Section 1.3.2 while ab initio simulations based on DFT have a limited scope in terms of the number of atoms simulated, they explicitly account for electronic polarization effects that govern the ion adsorption behaviour, and model charge transfer that can occur between the surface and the ions.
To overcome limitations on the number of particles in their simulations, Zhan et al. made use of an implicit solvent model, which replaces the water molecules with a background dielectric continuum. 143 Specifically, they use the reference interaction site model (RISM) that is capable of describing the electrolyte molarity using an effective screening approach; 180 this permits them to model the properties of the interface between graphene and a 1 M electrolyte solution using just a single adsorbed ion. Zhan et al. considered the different adsorption behaviour of Li + , Na + , K + and Cs + at a charge-neutral and negatively charged graphene electrodes. They found generally weaker interactions between the surface and the ions than Williams et al., 51 which was attributed to the electrolyte concentration in the RISM approach, compared with the pure water CPCM method, However an additional factor could be the inclusion of the point charge-quadrupole moment interactions that arise from the use of a hexagonal graphene flake in Williams et al.'s approach. 181 Zhan's work suggests that of the ions considered, only the K + and Cs + ions adsorb favourably on the surface, i.e. within the Helmholtz plane. Instead, they find that, from a thermodynamic perspective, the Li + and Na + ions prefer to remain in solvated in bulk, which is in keeping with our QMMD results. 40 However, it is observed that the minimum in the potential energy surface increases in distance from the graphene surface with the increase in ionic radius. Zhan et al. indicate that the negative adsorption energies of the larger ions, further away from the surface are not just due to the hydration free energy difference between the ions as asserted from classical simulations, but also due to a charge transfer mediated stabilization mechanism. They performed Bader charge density analysis 74 on the optimized electron density distribution in order to infer the fraction of charge transferred from the graphene to the adsorbed ion. This analysis revealed that even at the closest adsorption distances, there is negligible charge transfer between the smaller ions (Li + and Na + ) and the surface. However, the larger ions, with considerably more diffuse valence states that can strongly overlap with the p-electron cloud of the graphene, exhibit a significant amount of back electron transfer in the direction from the graphene to the cation. For example, the Cs + ion 3 Å (optimal computed adsorption height 3.5 Å) from the surface exhibits a computed partial charge of around 0.9 e, indicating a 10% fraction of electron transfer. This large fraction of charge transfer is thought to contribute to the energetic stabilization of the ion on the surface. These results underline the value of electronic structure theory and the necessity of multiscale simulations in the determination of interfacial properties.
A clear limitation of electronic structure theory based approaches is the computational cost associated with explicit solvation of the ions. Much of the available work describes ions which bind very strongly to the graphene electrode surface, both in the case that graphene is modelled as an infinite sheet or by a finite structure. Yet contrasting these results with Classical molecular dynamics, it is hard to ignore the fact that few classical models predict adsorption of the ions on the surface due the presence of explicit water. Ab initio molecular dynamics simulation of the graphene-electrolyte interface would address this issue, however obtaining sufficient time scales and ionic concentrations represents a severe challenge to this approach. With improved software algorithms, one avenue of exploration is to apply static DFT simulations directly on top of snapshots from Classical MD trajectories, 182 this has recently been applied to the graphene-ionic liquid interface by Da Silva et al. 183 and could form the basis for future investigations of the graphene-electrolyte interface.
Finally, based on models outlined in Section 1 and on our QMMD results 40 we have seen that the different adsorption behaviour of ions at the interface can underpin the EDL capacitance and device performance. Other MD simulations of the EDL capacitance have been conducted, aiming to reveal the role of specific ion effects and in particular the role those different cations have on the measured value, with somewhat conflicting results. Jiang et al. examined the EDL integral of capacitance Na-, K-, Cs-and RbCl salts by means of MD simulations with no polarization of the surface. 118 Contrary to our finding that the integral capacitance is constant, Jiang et al. report a decreasing C EDL with increasing ionic radius. Dočkal et al. also very recently computed the integral capacitance of the Chloride salts parameterized by Williams et al. 51,146 Their results again contradict previous simulations suggesting that the C EDL increases with ionic radius. Yang et al. also performed MD simulations of the same chloride salts, as well as mixtures of the different cations and were unable to detect any change in the computed capacitance. 184 One key difference between our work 40 185 This change to the dielectric constant close to the graphene has a direct impact on the computed potential drop and therefore can lead to the different trends observed in the simulations. In order to dispute or validate the results of these simulations, further characterization of the interface i.e. by techniques such as in-operando X-Ray absorption spectroscopy (XANES/EXAFS) or X-Ray reflectivity which provide information on the local structuring of the ions are required.

Effect of confinement on the electrolyte properties and its dynamics
Understanding the thermodynamics and dynamical properties of electrolytes under confinement is paramount for the design of carbon-based devices such as capacitors or membranes for water treatment. Graphene has been increasingly used as model material to build nanopores with controlled geometry and size to unravel the complex interplay between surface effects and bulk properties occurring at the nanoscale and how this affects important physical properties such relative permittivity, friction, surface tension and diffusion. [186][187][188][189][190] Carbon nanotubes (CNT) are probably the first example of carbon-based nanopores with well-defined sizes. The highly tuneable diameter of the tubes allowed controlled experiments to be performed that (in conjunction with simulations) indicated that narrow pores are characterized by extraordinary transport properties in terms of efficiencies and selectivities. [191][192][193] More recently a new technique pioneered by Radha et al. 194 opened up the possibility to fabricate atomically smooth two-dimensional graphite channels by removing selectively individual graphene sheets from graphite bulk crystals. The ability to build such highly size-controlled 2D capillaries has allowed direct measurements of the relative permittivity of water under confinement 195 and to characterize the water and electrolyte flow properties. 196 A direct comparison between the flow data obtained for the graphite capillaries and the CNTs is difficult as many of the interfacial properties of the latter seems to be highly dependent on the tube diameter (i.e. be curvature dependent). 197 Below we focus only on data obtained on the former.
Permeation experiments using both external electric fields or concentration gradients have shed light on the dynamics of the ions within 2D capillaries and contributed to develop better material-specific continuum theories. Since the early work 198 the permeation data showed a certain degree of ion specificity indicating that the ions' dynamics are not only driven by their bulk properties (i.e. hydration radius) but also by surface effects. Esfandiar et al. 198 measured the ionic conductivity of KCl solutions of different concentrations (between 1 and 10 À6 M) in graphene capillaries as small as 0.67 nm in height. They found that at high salt concentrations the ionic conductance was proportional to the electrolyte concentration and matched the KCl bulk conductivity, but for lower concentrations, 10 À4 M ca., the linear proportionality ceased to exist. Their graphene capillaries show a very low surface charge (r r 20 mC m À2 , order of magnitudes smaller compared to carbon nanotubes) and the dependence of the conductivity with the salt concentration was then attributed to OH À adsorption on the capillary walls. The authors also found that exchanging the K + ion with cations with larger hydration radius (larger than the capillary size) the electrolyte conductivity would reduce but would not be zero. They concluded that either the cations solvation shell would deform, or the ions would partially dehydrate. Finally performing drift-diffusion experiments at constant electrolyte concentration, they also noticed that, unlike K + ions, the Cl À ions have lower mobility under confinement than in bulk and also lower mobility than several of the cations investigated, despite having similar or smaller hydration radius. They explained this effect again in terms of the presence of OH À groups on the capillary walls that, interacting with the water in the Cl À solvation shell, would slow down the ions. The same group 199 manufactured graphene-slits of height of 0.34 nm and show that with such capillaries size only the proton could permeate.
Jung et al. 200 investigated the ionic conductivity of deionised water and NaCl solutions of 0.02 and 0.2 M concentrations in graphene and SiO 2 capillaries of sizes 3.6, 10 and 50 nm. 200 They found that graphene nanochannels show significantly higher ionic conductance than SiO 2 ones at all conditions investigated but the enhanced conductivity was more dramatic in the smaller channels and decreased with the salt concentration (see Fig. 8). The adsorption of hydroxide ions onto the surface of the graphene capillaries (in this case to explain the higher permeation of the ions in the graphene channels) and a change in slip length, are proposed to explain the difference in ionic mobility data between graphene and SiO 2 channels. The effect of the ionic concentration was explained in terms of Debye length which, decreasing with increasing concentration, was assumed to become smaller than the channel height for the smallest capillaries, thus increasing the electroosmotic ion transport and plug flow. The overall ion mobility increased with the channel sizes but the difference between the ionic conductivity measured in the graphene and silica nanochannels reduced with the channel size.
The experiments from Cheng et al. 201 clearly showed that conventional continuum theory cannot be used to predict the ion mobility for electrolytes under extreme confinement due to modification in the electrical double layer (EDL) structure. The authors fabricated graphene capillaries with interlayer distances between 0.8 nm and 5.4 nm and measured the ions permeation of KCl and K 2 SO 4 solutions of different concentrations as a function of the EDL thickness. The latter was tuned by changing the surface potential either by varying the externally applied potential or the degree of nanoconfinement. They found that both the potential sign and the degree of confinement affect the ion permeability which is ion-specific: applying a small negative voltage enhanced the ionic flux by up to 4 times compared to just a minor increase if a positive voltage of the same magnitude was used. Their experiments also showed that the ion permeation rates could be spontaneously changed, depending on the potential applied. The enhanced permeation was evident at strong confinement (capillaries smaller than 2 nm) and low electrolyte concentrations, but it was smaller for the larger channels (5.4 nm). The authors interpreted the results in terms of EDL atomic structure and the water arrangement in the first solvation shell of the ions. The fact that the potential sign affected the permeation indicated, the authors claimed, that the EDL structure differs depending on the nature of the counterions. This was confirmed by comparing the results across different cations and anions, which indicated a significant role in the permeation data of specific noncovalent interactions between the hydrated cations and anions in the EDL structure modulation. The authors then modified the Poisson-Nernst-Planck model to take into account that within small channels oscillating counter-and co-ion peaks in the ions' density profile should replace the ideal, exponential decay of the ion concentrations in the EDL predicted by the theory in bulk (see section 2.1). The authors showed that using a new analytical model that includes the ions' correlation effect, the trend of the calculated flux data with electrolyte concentrations and degree of confinement agree with their experiments.
This ions-dependence and flux enhancement, when a small external voltage is applied, as observed by Chen et al. has been confirmed by Mouterde et al. 202 who fabricated slit-like channels of both graphite and hexagonal boron nitride (hBN) of height 0.68 nm. The authors found again that Cl À anions have lower mobility than K + cations and that, in the graphitic channels, even with no external potential applied, the streaming mobility (i.e. the pressure driven component of the ionic current normalized by the channel length, pressure gradient and channel cross sectional channel) of potassium is higher than that reported for SiO 2 channels. When a small positive voltage was applied (75 mV), the potassium mobility was increased by 20 times. This gating effect was observed in both graphite and hexagonal boron nitride channels, but it was remarkably larger in the former. The authors attributed this high ionic mobility mainly to the fast transport of water (and hence hydrated ions, at molecular distances from the channel surfaces) observed in these types of channels and to the existence of an effective water-wall friction that depends on ion concentrations rather than to a change in the wall (graphite or hBN) capacitance. They then developed an 'Extended Poisson-Nernst-Planck' model that qualitatively reproduced most of their experimental observations including the large increase in streaming mobility under applied bias.
Due to the difficulty associated in building such nanoscopic channels and the possible contamination of the capillaries, most of the literature deals with molecular simulations work where it is much simpler to build perfectly aligned, rigid, graphene sheets sandwiched to form 2D symmetric channels of different heights, with the aim of disentangling the various competing effects affecting the permeability data. Molecular simulations have helped in clarifying important aspects of the ions mobility under confinement, for example that the ions rejection in small capillaries is dictated by the ion's hydration free energy rather than its van der Waals size, 203,204 that the preferential adsorption of some ions onto the graphite wall reduces the free energy of permeation into the capillaries 205 (see also Section 3.2), that water molecules when intercalated within the graphitic capillaries form discrete layers with lower dielectric constant 206 and that both Coulomb and Lennard-Jones interactions play a role on the ions' adsorption on the capillary walls depending on the channel size. 207 Molecular simulations have also been used to investigate the reduction in the dielectric constant of the liquid trapped inside nanocapillaries as its value affects the capacitance (see eqn (7)), EDL structure and ultimately flux data. 16,206,207 Experimentally Fumagalli et al. 23 measured the e r of water trapped in 1 nm graphite capillaries and found it to be as low as 2 (compared to 80 in bulk). However due to the challenges in carrying out the measurements and in the difficulty of ensuring no hydrocarbon contaminants are present in the channels, experimental values are scarce. Simulations have tried to shed light on the physical origin behind this reduction which happens also in more ''soft confinement'' than carbon slits. 22,208 Most of the literature has focused on water and generally ascribes the low e r calculated in the simulations to the molecular ordering and related reduced mobility of the confined water molecules (''ice-like'') particularly in the direction orthogonal to the channel walls, 16,22,[207][208][209] although a recent work by Olivieri et al. 185 proposed a new interpretation of the data and argued that the low dielectric constant value arises from the anisotropy of the long-ranged dipole correlation combined with an excluded-volume effect of the low-dielectric confining material.
Computer modelling allows also to design specific geometrical arrangements of the graphene sheets which might be difficult to achieve experimentally but nevertheless can provide guidelines for device design to maximize for example desalination processes or energy storage. 210,211 Of course, the availability of reliable molecular model parameters to simulate water and electrolytes under confinement is vital to obtain useful results. For example it has been shown that different water models have different phase diagram, 122,212,213 friction coefficient 214 and viscosity 215 within sub-nanometer carbon slits. The importance of reaching a consensus within the community on the most accurate model to be used for water under confinement and for the graphene/ electrolyte interactions (hampered so far also by the scarce number of experimental data) is proved by the plethora of (sometimes contradictory) results published in the literature including on the structure of the water. Moreover, as the model parameters employed in the simulations are different, it is difficult to paint a single picture of the behaviour of confined electrolyte. In general, molecular simulations indicated that for some degree of confinement and for some electrolytes, the ions mobility can be increased by confining the liquid. One of the early molecular dynamics study of NaCl electrolyte confined in infinite carbon slits 216 found that the diffusion coefficient of the ions increased by up to 20% in the largest slit (4 nm) while it was the same as bulk or slightly reduced for smaller interlayer distances (down to 1 nm). The same authors went on to investigate several other electrolytes of monovalent ions 217 (LiCl, NaCl, KCl, RbCl, CsCl, NaF, NaCl, NaBr, NaI) in the 4 nm slit and found that ion pairing either increased or decreased with confinement depending on the counterions and the anions without following a specific trend. Regarding the ions mobility the authors found that the chloride diffusion coefficient did not change substantially irrespectively to alkali ions but that potassium had significant rise in diffusivity relative to the other cations studied. They attributed this result to the lowering of water-water hydrogen bonds relative to the sodium chloride solution in the slit.
Chialvo and Cummings 140 simulated electrolytes of metal chloride (Li + , Ba 2+ , Y 3+ , Cl À ) in free-standing graphene slits with interlayer distances up to 3 nm. This set-up (i.e. open capillaries in equilibrium with a reservoir) had the advantage, compared to the infinite channel set-up, of removing any bias in terms of electroneutrality or density. The authors found that letting the ions move freely from the channel to the reservoir and vice versa show that as the graphene plate separation became smaller than the size of the ions first solvation shell, the ion would leave the capillary and move into the surrounding bulk electrolyte effectively charging the capillary. They also found that the slit-pore would de-wet as the interplate separation is reduced below 0.9 nm ca.
Sala et al. 218 focused on aqueous solutions of sodium halide (NaF, NaCl and NaI) of different molalities (1, 2, 4 mol kg À1 ) confined inside a graphene channel of size 3.1 nm. Unlike the previous studies the authors employed a polarizable model for both water and ions but not for the graphene. From a structural point of view, the effect of using a polarizable force field made a difference only for the position of anions (which are fully hydrated in the channel or closer to the capillary walls when modelled with non-polarizable or polarizable force field respectively). Both force fields predicted higher permittivity at the interface compared to the bulk value (in contrast with the standard theoretical and simulation results 144,219 ) and similar diffusion coefficients. In terms of ions relative diffusion, the effect of confinement was found to be minor for all the ions except for fluoride whose dynamics slow down.
Kong et al. 220 investigated the effect that temperature (283, 298, 313 and 333 K) and surface charge have on the ions dynamics simulating a NaCl solution confined in a graphene slit with heights of 0.7, 1.2, 1.7 nm, respectively, in contact with two 2 M electrolyte reservoirs. In this case the constituent atoms of the graphene have a small partial charge associated (the same for all carbons) in order to reproduce a surface charge of |5| mC cm À2 . As in the set-up of Chialvo, 140 the capillaries are in contact with a solution reservoir and therefore their electrolyte concentration may vary. The results show that ions confined in neutral nanochannels diffused faster than in charged ones (this was attributed to the interaction between the charged wall and the counterions) and the electrolyte conductivity in charged channels decreased with decreasing channel height. As expected, an increase in temperature led to enhanced thermal motion of ions, but the overall effect of the temperature of the diffusion was negligible.
Having the capillaries in contact with the electrolyte reservoirs allowed Kalluri et al. 221 to investigate the charging of graphitic slit pores of different sizes as a function of electrolyte concentration and surface charge. The pores of size 0.9, 1.2, and 1.6 nm and charge densities ranging from 0 (neutral pore), 20, 30, and 40 mC cm À2 were soaked in aqueous solutions of NaCl at 1.5 and 1.6 M. The authors found that the pores were fully wetted by the electrolyte solutions only when they were charged and that at the maximum graphene surface charge density considered, the ionic concentration within the pores could become B10 times as high as that outside the pores. Such high surface charge led to ions condensation on the slit walls sometimes in multiple layers. Within this set up and FF choice, capillaries of width 1.2 nm seem to be those with the highest charge accumulation.
The effect on the channel electrification on the capacitance and EDL structure has been investigated by Feng et al. 222 who modelled graphene slits with widths ranging from 0.94 to 1.47 nm and surface charge of À5.5 mC cm À2 . Only potassium cations were dissolved into water to neutralize the simulation box and no counterions were included. The MD simulation results showed that the cations concentration profile followed the classical EDL theory (i.e. adsorbed at the wall) only for the smallest slit width (0.94 nm) while for larger channels the ions tended to accumulate in the middle of the capillary. In the smallest channels however, they found that the adsorption of the cations to the capillary walls led to a depletion of ions in the middle of the channel. The authors attributed this result to both effects accounted for in the original EDL theory and others arising from the discrete nature of the electrolyte and therefore not included in the EDL. The former include the long-range electrostatic ion-ion repulsion (driving ions toward the two slit walls) and entropic effects (driving ions and water away from the wall); the latter are the van der Waals ion-slit wall attractions, the hydration free energy of the ions (driving ions into the middle of the channels to be fully hydrated) and the interactions between an ion's hydration water molecules and their surrounding water molecules. The ions/wall attraction due to surface polarization was excluded in this analysis as the molecular model employed did not include this effect which however exists (see Section 3.2 above). The authors then proposed to calculate the capacitance of nanopores approximating the slit capacitor as two plate capacitors in parallel and calculated the differential Helmholtz double layer capacitance of each wall using eqn (5) (rather than using the ions distribution profile obtained via the MD simulations and eqn (18)). Despite this approximation assuming a specific geometry of the pore, they then showed that the parallel plate capacitor approximation fitted the experimental data well for the capacitance of activated carbon.
More speculative work, which is harder to compare or validate with the experiments due to the greater mismatch between the idealized model and the actual device includes the calculations of the graphene edge's capacitance. MD simulations have been performed on aqueous electrolytes in contact with vertically-oriented graphene flakes. 223,224 The results of these simulations need to be taken with caution though, as the H atoms of the graphene edges are much more reactive than the carbon atoms of the basal plane and therefore the assumption made in the model that the edges are not functionalized makes the comparison with the experimental data difficult. Anyhow, the MD results indicated that the specific capacitance, C diff , of the edges can achieve a 2-fold increase compared to the basal one when the interlayer spacing is B0.5 nm. An analysis of the EDL structure showed that the calculated capacitance arises from the presence of the water in the EDL rather than the adsorption of the ions. Experimentally the edge capacitance is also larger than the basal one but the increase is much higher than the simulation results indicated, probably due to the functionalization of the graphene edges which is not taken into account in the model. 92,225

Conclusions and outlook
Despite the intense research done on the topic, a molecularlevel understanding of the electrochemical properties of graphene (including few layer graphene) and aqueous electrolyte is still elusive. The possibility of using experimental methods in conjunction with simulations, represents a powerful tool that can be exploited further as long as experiments and simulations look at the same system. On the simulation side, considerations of the electrochemical stability of the electrolyte solutions at high applied voltage should be given particular attention and therefore simulations with high surface charge should be avoided. If electron transfer is accurately included, however, then simulations could offer valuable insight into the range of electrochemical stability of graphene vs. graphite in different electrolyte solutions. Simulations have also quite clearly indicated now that including polarization in models of the electrolyte solution has a minor effect on the structure of the interface but including the surface polarizability is important to simulate the ion's adsorption and wetting. The routine use of defect-free pristine graphene in the simulation should also be reconsidered or at least tested as even the most controlled experimental set-ups present surface defects and, most probably, contaminations including the presence of various functional groups. The latter observation is even more pertinent when simulations are carried out on graphene edges whose capping hydrogen atoms are particularly reactive. The surface defects, in particular for graphene, could also have particular relevance in the quantum capacitance calculations as the latter is intimately related to the DOS of the material. As discussed earlier in the review, a range of possible structural defects can exist, and the experimental complication is the difficulty in ensuring that only one class of defect is introduced into the graphene structure. This is, again, where interplay between experiment and simulation becomes particularly useful, because the latter can indicate which defects would be most likely to produce the largest perturbation of the sample's DOS and hence capacitance. A further effect which has seen relatively little exploration from the experimental side, is the extent to which the underlying substrate influences the properties (e.g. capacitance) of the graphene/electrolyte interface. Monolayer graphene samples are required to maximise the strength of this phenomenon.
Furthermore, for a quantitative comparison between simulation and electrochemical experimental data, several points need still to be clarified. For instance, it is still unclear what value of the local permittivity should be used for the capacitance calculation (eqn (7), (16) and (17)) considering both the higher ionic concentration in the EDL compared to bulk and the effect of the proximity of the surface 185 or if interfacial crystallization occurs in the EDL layer 144 opening up the possibility of a completely different interpretation of the experimental results. Simulations are particularly important for shedding light on the effect that confinement has on the electrochemical properties. As said, much of the work has been focused on clarifying the atomistic structure of the confined EDL but information about the electrolyte physical properties beyond the EDL are also vital. These parameters are important because they inform continuum equations to predict ion and water diffusivity and channels flux properties. Among them wetting phenomena (relevant also in the bulk) occurring in nano-confined systems 189 although likely to be relevant, have not been heavily explored.
Finally, this review has focused at low-medium concentrated electrolytes, but a new type of highly concentrated electrolyte (known as water-in-salt) 226 has been emerging as environmentallyfriendly, water-based electrolyte for use in energy storage devices. These seem to represent a new frontier in the field and much of their electrochemical properties are still to be explored in detail. For example, the general trend observed is that these highly soluble salts greatly extend the stable potential window at the graphite (carbon)/aqueous electrolyte interface, to values way in excess the thermodynamic 1.23 V discussed in section 1.5 above. The factors underlying this stabilisation are not entirely clear, with a reduced thermodynamic activity of water (associated with the loss of ''free'' water molecules) being an obvious explanation, however some works claim that at least part of the stabilisation is due to the decomposition of some of the ionic species to form a solid-electrolyte interface layer, which insulates the water from the electrode and thus suppresses the hydrogen evolution process. 226,227 Understanding how graphene interacts with such highly concentrated electrolytes will be extremely interesting, both from the fundamental and ''technological'' perspectives.

List of symbols and abbreviations
A Surface area of the electrode c Electronic charge carrier density in graphite C Capacitance of the parallel plate capacitor C diff Specific differential capacitance C d,eff Effective dense layer capacitance

Conflicts of interest
There are no conflicts to declare. 760907 and the Leverhulme Trust (grant RPG-2019-250) for funding.