Do specific ion effects determine performance of aqueous graphene-based supercapacitors? Perspectives from multiscale QMMD simulations

Whether or not specific ion effects determine the charge storage properties of aqueous graphene and graphite based supercapcitors remains a highly debated topic. In this work we present a multiscale quantum mechanics – classical molecular dynamics (QMMD) investigation of aqueous monoand divalent salt electrolytes in contact with fully polarizable charged graphene sheets. By combining the electrochemical double layer (EDL) and quantum capacitance we observe an increasing electrode specific capacitance with cationic radii and charge. Counterintuitively, we determine that a switch in the cation adsorption mechanism from inner to outer Helmholtz layers leads to negligible changes to the EDL capacitance. However, the ability of ions (such as K) with a relatively low hydration free energy to penetrate the inner Helmholtz plane and adsorb directly on the electrode surface is found to slow their diffusion dynamics. Our results show that surface effects underpin the behaviour of cations at the interface and add a vital new perspective on trends in ion mobilities seen under confinement. Aqueous carbon-based supercapacitors represent a chemically stable, low-cost and non-toxic solution to low charge storage capacities and slow charge/discharge rates in conventional capacitors and batteries. Graphene, graphite and carbon nanotube electrodes each offer comparatively high surface areas that increases the number of ionic species that can be stored electrostatically within the electrochemical double layer. The reversible, non-Faradic physisorption of ions at the electrode surfaces mediates a rapid release of stored charge and therefore drives the function the supercapacitor. The challenge that remains is to identify and fine tune exact electrode-electrolyte compositions that simultaneously maximise the density of surface ion adsorption and minimizes the charge/discharge rate. Critical to the continued development of aqueous carbon-based supercapacitors is an understanding whether or not ion specificity plays a role in the amount of charge stored at the interface, and if so: to what extent. Currently, experimental measurements in this area are conflicted, with mixed reports on the behaviour of different ions. The capacitance obtained from graphene films (synthesised from graphene oxide) in contact with LiCl, NaCl and KCl supports the idea that water as a polar solvent inhibits any specific ion effects. On the contrary, in a separate investigation of the capacitance of graphite electrodes, a survey of group 1 chloride salts found that while specific ion effects could not be detected for the edge-oriented graphite, a clear increase in the capacitance was measured with increasing ionic radii. These specific ions effects were linked to the decreasing free energies of hydration as the cations get larger. Atomistic simulations based on either fully classical models or electronic structure theory provide an atom-scale resolution of the structuring and staticand dynamical properties of the electrode-electrolyte interface. These simulations increase in computational complexity from classical molecular dynamics (MD) that uses force fields to model the interactions between atoms to first principles molecular dynamics (FPMD). Specifically in the investigation of the charge storage mechanism in aqueous supercapacitors, the non-bonded electrostatic interactions brought about by ion (or water) induced polarization of the electrode surface can be critical in reproducing experimental observations. These interactions are explicitly included in FPMD, which enables investigation of the properties of many different aqueous electrode interfaces. However, large scale models that capture the entropic effects, known to be important for interfacial properties, require simulations of several thousands of atoms for tens of nanoseconds, which for FPMD methods are prohibitively expensive. This means that a balanced accuracy-computational viability trade-off must be established. Parameterization of the MD non-bonded potential to account for polarized ion-surface interactions facilitates simulations on the necessary time and length scales, with the caveat that the cumbersome parameterization must be carried out for each and every ionic species, solvent model, surface morphology and surface charge density. In addition static nature of polarization contained in LennardJones potential overpredicts binding energy by more than 10 kJ mol. Alternatively, a description of the dynamical surface polarization can be introduced into the MD force field by tethering a dummy charge to the atom via a harmonic spring. Similarly, this type of polarizable force field requires a careful parameterization for scenarios where the surface electrode is charged, but also limits the polarizability of the surface to localized spatial regions in the vicinity of each atom and therefore potentially misses long-ranged polarization that results from redistribution of electronic charge density. For strictly metallic surfaces, where the potential in the electrode is fixed, the constant potential method modulates the electrode atom partial charges in response to the specific interfacial electrolyte configuration. The validity of this approach for semimetallic electrode surfaces was recently called into question. On the other hand, quantum mechanical (QM) simulations of the ion-surface interactions can be carried out in the presence of an implicit solvent model. Removing explicit water molecules drastically reduces the computational cost of the QM simulations and in addition, the treatment of ion and surface electron densities also introduces charge-transfer (CT) interactions, that can contribute to the stabilisation of ions at the surface, to the model. Yet, this comes at the expense of a lack of both the interface dynamics and the explicit ion-hydration shell that contribute to specific adsorption behaviour. Leveraging the explicit polarization of QM approaches and the rapid simulation of large length and time scales in the classical MD regimes, we previously introduced an iterative QMMD scheme that embeds the dynamical electronic structure of the electrode within a classical force field. This method retains a description of the long-ranged polarization of the surface electrode irrespective of its metallicity, while allowing for nanometre and nanosecond simulations. However, fully classical treatment of the electrolyte phase necessarily means that CT interactions that occur on the Angstrom scale are not included. In this work we deploy our QMMD framework to model the capacitance of graphene electrodes as depicted in Figure 1 at three fixed surface charge densities: the neutral case and charged negatively and positively with ±0.061 Cm. We consider 1.0 M solutions of LiCl, NaCl, KCl, MgCl2 and CaCl2 recently parameterized for the TIP4P/2005 water model using scaled ionic charges. Further simulation details are provided in the Supporting Information. We uncover a switch in the monovalent cationic-adsorption mechanism as the ionic radius increases, which helps to resolve outstanding controversies in specific ion effects at the graphene interface. Figure 1: (Top) Representative simulation box for the graphene-aqueous electrolyte interfaces considered in this work. The two-dimensional graphene electrode (cyan), approximately 3 × 3 nm in dimension, is placed in contact with 8 nm slabs of electrolyte (cations: purple, anions: yellow, O: red, H: white) with a further 8 nm vacuum buffer to prevent interactions between periodic images. (Bottom), Plot of the bulk normalized densities along the surface normal for each of the different electrolytes considered for a neutral (left) and negative (centre) and positively (right) charged electrode. Vertical dashed lines denote the positions of the first cation adsorption peaks. We first examine the structure of the various graphene-electrolyte interfaces at 1M concentration, including the effects of the quantum mechanically polarized surface. The number density profile along the axis normal to the surface describes the coarse structuring of the electrochemical double layer at the interface. To extract information from the density profile we introduce the bulk-normalized number density, ?̅?(z), which is scaled so that that density for each given species in the bulk region is unity. The ?̅?(z) are plotted in Figure 1. Across all the different electrolyte solutions, the water (blue, Figure 1) and Cl anion (yellow, Figure 1) density profiles do not change appreciably; the respective first adsorption peaks are located at 0.33 and 0.65 nm. The interfacial water adopts the classic structure, having two separate peaks; the placement of these peaks is consistent with previous investigations of graphene and the TIP4P-2005 water model. We note that the Cl anion is repelled from the neutral surface, a finding consistent with conclusions of Misra et al that strongly hydrated kosmotrophic anions are generally repelled by graphene. In the case of the monovalent cations, the ?̅?(z) of the smaller Li and Na indicate that these ions reside in the outer Helmholtz plane, likely due to their respective high free hydration energies and kosmotrophic behaviour preventing partial dehydration of the ion first hydration sphere, this is discussed in greater detail below. In the case of Li, the first adsorption peak is found at 0.41 nm and the second at 0.72 nm, whereas Na has a minor adsorption peak at 0.37 nm, a major adsorption peak at 0.51 nm and a second adsorption peak at 0.75 nm. Each of these peaks is further away from the electrode than first water layer, which is consistent with our previous finding for Na using the QMMD approach, which incidentally used different classical parameters. Like the Na ion, the larger K cation (purple, Figure 1) also has minor and major adsorption peaks, these are located at 0.29 nm and 0.58 nm, with a second adsorption peak at 0.79 nm. It is noteworthy that the minor peak is closer to the graphene surface than the first water peak indicating that the K ion can penetrate the inner Helmholtz plane. This is at variance with the finding of Pykal et al, who investigated aqueous KFand KI-graphene interfaces using a polarisable force field model, and always observed K in the outer Helmholtz or diffuse layer. In their work it was deemed that the K ion adsorption behaviour was linked to the adsorption behaviour of the anion, with iodide able to displace K at the surface and fluoride able to draw the K cation into the bulk. It is therefore conceivable that the position of the Cl ion in the outer Helmholtz layer facilitates adsorption of K at the interface. However this could also be linked to the ability of our model to delocalize surface charge density to the adsorption site, an effect not captured with the Drude oscillator polarizible force field. We highlight that the large K cation is generally considered to be chaotrophic and therefore the adsorption of K in the inner Helmholtz layer is in further agreement with similar conclusions drawn for anions, our findings suggest that chaotrophic cations are also attracted to the neutral graphene surface. While there is diversity in the profiles of the monovalent cations, the two divalent cations, we consider Mg (red, Figure 1) and Ca (magenta, Figure 1) have highly similar density profiles: the first and second adsorption peaks are located at 0.47 & 0.73 nm and 0.49 & 0.78 nm respectively, in the outer Helmholtz and diffuse layers. In addition to the examination of the behaviour of the system setup with a charge neutral electrode, we also computed the ?̅?(z) for 1.0 M electrolyte concentrations in contact with electrodes charged positively and negatively, these are also presented in Figure 1. In order to maintain overall electroneutrality, the electrode surface charge density has been set to ±0.061 Cm and the number of ions in solution modulated to balance the electrode charge. Our results show that upon charging, only the K ion enters the inner Helmholtz layer and directly adsorbs on the graphene surface. When the graphene electrode is charged negatively the positions of main peaks in the ?̅?(z) profile relative to the graphene electrode remain the same. However, the respective peak intensities, specifically those at the interface, vary drastically. At 1.0 M concentration, the first and second adsorption peaks for Li are again found at 0.4 and 0.7 nm from the graphene surface. As with the neutral electrode, these are still in the outer Helmholtz plane since the first water peak is at 0.3 nm, but now with respective intensities 5 and 1.5 times larger than the bulk density. The same effect is seen when we consider Na, the first and second adsorption peaks are observed 0.5 and 0.7 nm respectively, the peak intensities are 3 and 1.5 times larger than the bulk density. On a technical note, the fact that the cation peaks do not change position upon negative charging of the electrode indicates that the in the trade-off between the repulsive part of the Lennard-Jones interaction potential and the attractive surface polarization-augmented Coulomb potential, the Lennard-Jones interaction dominates. This can be explained from the results of Misra et al, who identified that interference between water and ion electric fields effectively screens up to 85 % of the surface ion attractions. The small difference in the positions of the first adsorption peak of the Li and Na (Na 0.1 nm further from the surface) is likely due to the larger ionic radius and solvation sphere of the Na ion. The result is a slightly weaker attractive Coulomb interaction with the charged surface and therefore reduced adsorption intensity; averaged over 15 ns of QMMD simulation the short-ranged Coulomb interaction between the anode and ion is -28.4 kJ/mol in the case of Li and -20.9 kJ/mol for Na. As shown in Figure 1, for K the intensity of the first peak at 0.29 nm, which is within the inner Helmholtz plane, is six times greater than the bulk concentration indicating a significant build-up of ions in direct contact with the negatively charged electrode. Consequently, we observe that the K ion has a strong Coulombic attraction to the polarizable surface -74.6 kJ/mol, approximately three times greater than the other monovalent ions in the outer Helmholtz layer. Turning to the divalent ions, the first and second adsorption peaks are again found at 0.47 & 0.73 nm for Mg, and 0.50 & 0.76 nm for Ca. The intensities of the first adsorption peaks are five and four times the bulk ion density respectively, which further indicates that there is little to discriminate between these two ions. The Coulombic attraction between the divalent ions and the electrode is expectedly larger than the monovalent ions Li and Na, Mg -40.5 kJ mol Ca -30.4 kJ mol, this is due to their similar proximity to the surface and larger ionic charge. Charging the electrode positive has the same effect on all the systems considered. The cation first adsorption peaks closest to the electrode are reduced to zero, this is due to positively charged ions being repelled. Instead, we observe the outer Helmholtz layer adsorption of the Cl anions at approximately 0.4 nm from the surface with an intensity two times the bulk density. In order to determine which properties are the driving force for the K ion to penetrate the inner Helmholtz layer, in contrast to the other cations considered, we performed a series of simulations with a negative charging of the electrode. In each simulation one cation was frozen either in the bulk solution or close to the surface; the surface absorption height was chosen according to the first adsorption peak in the density profile as reported in Table 1. Figure 2 reports the total number of coordinated water molecules within a given radius in bulk (a) and adsorbed on the surface (b). Within the first coordination sphere, which is 0.35 nm, the difference between the bulk and surface is minimal (< 0.1 water molecules) for all ions except K. The latter has on average 1.5 less coordinated water when it is adsorbed at the surface. These results are highlighted in Figures 2c and 2d, which show the coordination of a Na and K ion respectively at the charged surface. For Na, the closest 5 water molecules make up the first coordination shell, of which two clearly reside in between the ion and the surface. In contrast the K ion is directly adsorbed on the graphene surface with no intermediate water molecule. We note that even in the bulk K does not have a stable first solvation shell; unlike each of the other ions the K radial distribution function has a broad peak with an elongated tail which signifies that the final molecule in the coordination sphere is weakly bound at a moderately large distance from the ion. Consequently, there is a comparatively low energy requirement to remove such a molecule from the solvation shell with respect to the other ions investigated that permits the dehydration we observe at the surface. Behaviour which is consistent with results from investigations of K ion migration through biological channels. Electrolyte dBulk (nm) dSurface (nm) LiCl 3.991 0.414 NaCl 3.960 0.393 KCl 4.328 0.280 MgCl2 4.256 0.474 CaCl2 4.207 0.487 Table 1: Vertical distance between the electrode and the frozen cation in the bulk electrode and at the interface. Surface ions have been fixed at the vertical coordinate corresponding to the first peak in the density profile. Figure 2: Computed radial distribution functions and integrated radial distribution functions of frozen cations in the bulk solution (a) and at the surface (b) for the various 1M concentration electrolytes considered in this work. Representative configurations of the Na (c) and K (d) ions frozen at the interface. Direct exposure of the K ion to the surface leads to an induced polarization of the C atoms close to the ion resulting in strongly negative attractive Coulomb interactions. In Figure 3a, for a single example configuration, we plot the surface charges in two dimensions and mark the positions of K ions directly adsorbed to the graphene with crosses circled blue. As a guide to the eye, the surface charges are reported normalized to the homogenous fixed charge −3.4 336 ⁄ e (where −3.4e is the total charge of the graphene sheet made of 336 carbon atoms). In this example the surface Mulliken charges are significantly more negative in the surface region local to the ion. To understand the extent to which this is a general property of the directly adsorbed ions, we also plot in Figure 3b all surface charges as a function of the distance from each adsorbed K ion across the 15 ns trajectory. As seen by the linearly increasing trend, up to 0.35 nm the polarization of the C atoms is very strongly influenced by the ion, this corresponds to nearest and immediately bonded atoms in the surface. Beyond 0.35 nm the computed Mulliken charges average out to the homogenous fixed charge value through interactions with other ions or water molecule dipoles. Our results demonstrate the importance of the local ion induced polarization in the stabilization of the direct adsorption configuration and the clear trend and distance dependence opens up for the possibility of an analytical description of the polarization energy to further expedite QM simulations. Figure 3: (a) Representative plot of the computed Mulliken charges on the graphene sheet in contact with a 1M solution of KCl and charged with 3.4 e. Circled X’s mark the coordinates of K ions directly adsorbed on the surface. (b) A plot of the computed Mulliken charge as a function of the distance from each K ion adsorbed on the surface over an entire 10 ns QMMD trajectory. For illustrative purposes, all computed charges are normalized to the homogenously averaged surface charge -3.4/336 e. We turn now to the electrochemical double layer integral capacitance, which is obtained from the system number density profile along the surface normal, n(z) . First, the charge density, ρ(z) , is calculated from the atom-resolved number densities multiplied by their partial charges from the classical force field. This includes all charged atoms (appropriately scaled), ions and electrodes. Poisson’s equation relates the curvature of the electrostatic potential, Φ, to the charge density, which along one dimension reads,

Aqueous carbon-based supercapacitors represent a chemically stable, low-cost and non-toxic solution to low charge storage capacities and slow charge/discharge rates in conventional capacitors and batteries. Graphene, [1][2][3][4][5][6] graphite and carbon nanotube 7,8 electrodes each offer comparatively high surface areas that increases the number of ionic species that can be stored electrostatically within the electrochemical double layer. The reversible, non-Faradic physisorption of ions at the electrode surfaces mediates a rapid release of stored charge and therefore drives the function the supercapacitor. The challenge that remains is to identify and fine tune exact electrode-electrolyte compositions that simultaneously maximise the density of surface ion adsorption and minimizes the charge/discharge rate.
Critical to the continued development of aqueous carbon-based supercapacitors is an understanding whether or not ion specificity plays a role in the amount of charge stored at the interface, and if so: to what extent. Currently, experimental measurements in this area are conflicted, with mixed reports on the behaviour of different ions. The capacitance obtained from graphene films (synthesised from graphene oxide) in contact with LiCl, NaCl and KCl supports the idea that water as a polar solvent inhibits any specific ion effects. 9 On the contrary, in a separate investigation of the capacitance of graphite electrodes, 10 a survey of group 1 chloride salts found that while specific ion effects could not be detected for the edge-oriented graphite, a clear increase in the capacitance was measured with increasing ionic radii. These specific ions effects were linked to the decreasing free energies of hydration as the cations get larger. 11 Atomistic simulations based on either fully classical models or electronic structure theory provide an atom-scale resolution of the structuring and static-and dynamical properties of the electrode-electrolyte interface. These simulations increase in computational complexity from classical molecular dynamics (MD) that uses force fields to model the interactions between atoms to first principles molecular dynamics (FPMD). Specifically in the investigation of the charge storage mechanism in aqueous supercapacitors, the non-bonded electrostatic interactions brought about by ion (or water) induced polarization of the electrode surface can be critical in reproducing experimental observations. [12][13][14] These interactions are explicitly included in FPMD, which enables investigation of the properties of many different aqueous electrode interfaces. [15][16][17][18] However, large scale models that capture the entropic effects, known to be important for interfacial properties, require simulations of several thousands of atoms for tens of nanoseconds, which for FPMD methods are prohibitively expensive. This means that a balanced accuracy-computational viability trade-off must be established.
Parameterization of the MD non-bonded potential to account for polarized ion-surface interactions facilitates simulations on the necessary time and length scales, 19 with the caveat that the cumbersome parameterization must be carried out for each and every ionic species, solvent model, surface morphology and surface charge density. In addition static nature of polarization contained in Lennard-Jones potential overpredicts binding energy by more than 10 kJ mol -1 . 14 Alternatively, a description of the dynamical surface polarization can be introduced into the MD force field by tethering a dummy charge to the atom via a harmonic spring. 20 Similarly, this type of polarizable force field requires a careful parameterization for scenarios where the surface electrode is charged, but also limits the polarizability of the surface to localized spatial regions in the vicinity of each atom and therefore potentially misses long-ranged polarization that results from redistribution of electronic charge density. For strictly metallic surfaces, where the potential in the electrode is fixed, the constant potential method modulates the electrode atom partial charges in response to the specific interfacial electrolyte configuration. [21][22][23][24] The validity of this approach for semimetallic electrode surfaces was recently called into question. 13 On the other hand, quantum mechanical (QM) simulations of the ion-surface interactions can be carried out in the presence of an implicit solvent model. [25][26][27] Removing explicit water molecules drastically reduces the computational cost of the QM simulations and in addition, the treatment of ion and surface electron densities also introduces charge-transfer (CT) interactions, that can contribute to the stabilisation of ions at the surface, 28 to the model. Yet, this comes at the expense of a lack of both the interface dynamics and the explicit ion-hydration shell that contribute to specific adsorption behaviour.
Leveraging the explicit polarization of QM approaches and the rapid simulation of large length and time scales in the classical MD regimes, we previously introduced an iterative QMMD scheme that embeds the dynamical electronic structure of the electrode within a classical force field. 13 This method retains a description of the long-ranged polarization of the surface electrode irrespective of its metallicity, while allowing for nanometre and nanosecond simulations. However, fully classical treatment of the electrolyte phase necessarily means that CT interactions that occur on the Angstrom scale are not included.
In this work we deploy our QMMD framework to model the capacitance of graphene electrodes as depicted in Figure 1 at three fixed surface charge densities: the neutral case and charged negatively and positively with ±0.061 Cm -2 . We consider 1.0 M solutions of LiCl, NaCl, KCl, MgCl2 and CaCl2 recently parameterized for the TIP4P/2005 water model using scaled ionic charges. 29 Further simulation details are provided in the Supporting Information. We uncover a switch in the monovalent cationic-adsorption mechanism as the ionic radius increases, which helps to resolve outstanding controversies in specific ion effects at the graphene interface. We first examine the structure of the various graphene-electrolyte interfaces at 1M concentration, including the effects of the quantum mechanically polarized surface. The number density profile along the axis normal to the surface describes the coarse structuring of the electrochemical double layer at the interface. To extract information from the density profile we introduce the bulk-normalized number density, ̅( ), which is scaled so that that density for each given species in the bulk region is unity. The ̅( ) are plotted in Figure 1.
Across all the different electrolyte solutions, the water (blue, Figure 1) and Cl anion (yellow, Figure 1) density profiles do not change appreciably; the respective first adsorption peaks are located at 0.33 and 0.65 nm. The interfacial water adopts the classic structure, having two separate peaks; the placement of these peaks is consistent with previous investigations of graphene and the TIP4P-2005 water model. 30,31 We note that the Cl anion is repelled from the neutral surface, a finding consistent with conclusions of Misra et al that strongly hydrated kosmotrophic anions are generally repelled by graphene. 14 In the case of the monovalent cations, the ̅( ) of the smaller Li and Na indicate that these ions reside in the outer Helmholtz plane, likely due to their respective high free hydration energies 32 and kosmotrophic behaviour preventing partial dehydration of the ion first hydration sphere, this is discussed in greater detail below. In the case of Li, the first adsorption peak is found at 0.41 nm and the second at 0.72 nm, whereas Na has a minor adsorption peak at 0.37 nm, a major adsorption peak at 0.51 nm and a second adsorption peak at 0.75 nm. Each of these peaks is further away from the electrode than first water layer, which is consistent with our previous finding for Na using the QMMD approach, 13 which incidentally used different classical parameters. [33][34][35] Like the Na ion, the larger K cation (purple, Figure 1) also has minor and major adsorption peaks, these are located at 0.29 nm and 0.58 nm, with a second adsorption peak at 0.79 nm. It is noteworthy that the minor peak is closer to the graphene surface than the first water peak indicating that the K ion can penetrate the inner Helmholtz plane. This is at variance with the finding of Pykal et al, who investigated aqueous KF-and KI-graphene interfaces using a polarisable force field model, and always observed K in the outer Helmholtz or diffuse layer. 36 In their work it was deemed that the K ion adsorption behaviour was linked to the adsorption behaviour of the anion, with iodide able to displace K at the surface and fluoride able to draw the K cation into the bulk. It is therefore conceivable that the position of the Cl ion in the outer Helmholtz layer facilitates adsorption of K at the interface. However this could also be linked to the ability of our model to delocalize surface charge density to the adsorption site, an effect not captured with the Drude oscillator polarizible force field. 36 We highlight that the large K cation is generally considered to be chaotrophic and therefore the adsorption of K in the inner Helmholtz layer is in further agreement with similar conclusions drawn for anions, 14 our findings suggest that chaotrophic cations are also attracted to the neutral graphene surface. While there is diversity in the profiles of the monovalent cations, the two divalent cations, we consider Mg (red, Figure 1) and Ca (magenta, Figure 1) have highly similar density profiles: the first and second adsorption peaks are located at 0.47 & 0.73 nm and 0.49 & 0.78 nm respectively, in the outer Helmholtz and diffuse layers.
In addition to the examination of the behaviour of the system setup with a charge neutral electrode, we also computed the ̅( ) for 1.0 M electrolyte concentrations in contact with electrodes charged positively and negatively, these are also presented in Figure 1. In order to maintain overall electroneutrality, the electrode surface charge density has been set to ±0.061 Cm -2 and the number of ions in solution modulated to balance the electrode charge.
Our results show that upon charging, only the K ion enters the inner Helmholtz layer and directly adsorbs on the graphene surface. When the graphene electrode is charged negatively the positions of main peaks in the ̅( ) profile relative to the graphene electrode remain the same. However, the respective peak intensities, specifically those at the interface, vary drastically. At 1.0 M concentration, the first and second adsorption peaks for Li are again found at 0.4 and 0.7 nm from the graphene surface. As with the neutral electrode, these are still in the outer Helmholtz plane since the first water peak is at 0.3 nm, but now with respective intensities 5 and 1.5 times larger than the bulk density. The same effect is seen when we consider Na, the first and second adsorption peaks are observed 0.5 and 0.7 nm respectively, the peak intensities are 3 and 1.5 times larger than the bulk density. On a technical note, the fact that the cation peaks do not change position upon negative charging of the electrode indicates that the in the trade-off between the repulsive part of the Lennard-Jones interaction potential and the attractive surface polarization-augmented Coulomb potential, the Lennard-Jones interaction dominates. This can be explained from the results of Misra et al, who identified that interference between water and ion electric fields effectively screens up to 85 % of the surface ion attractions. 14 The small difference in the positions of the first adsorption peak of the Li and Na (Na 0.1 nm further from the surface) is likely due to the larger ionic radius and solvation sphere of the Na ion. The result is a slightly weaker attractive Coulomb interaction with the charged surface and therefore reduced adsorption intensity; averaged over 15 ns of QMMD simulation the short-ranged Coulomb interaction between the anode and ion is -28.4 kJ/mol in the case of Li and -20.9 kJ/mol for Na. As shown in Figure 1, for K the intensity of the first peak at 0.29 nm, which is within the inner Helmholtz plane, is six times greater than the bulk concentration indicating a significant build-up of ions in direct contact with the negatively charged electrode. Consequently, we observe that the K ion has a strong Coulombic attraction to the polarizable surface -74.6 kJ/mol, approximately three times greater than the other monovalent ions in the outer Helmholtz layer.
Turning to the divalent ions, the first and second adsorption peaks are again found at 0.47 & 0.73 nm for Mg, and 0.50 & 0.76 nm for Ca. The intensities of the first adsorption peaks are five and four times the bulk ion density respectively, which further indicates that there is little to discriminate between these two ions. The Coulombic attraction between the divalent ions and the electrode is expectedly larger than the monovalent ions Li and Na, Mg -40.5 kJ mol -1 Ca -30.4 kJ mol -1 , this is due to their similar proximity to the surface and larger ionic charge.
Charging the electrode positive has the same effect on all the systems considered. The cation first adsorption peaks closest to the electrode are reduced to zero, this is due to positively charged ions being repelled. Instead, we observe the outer Helmholtz layer adsorption of the Cl anions at approximately 0.4 nm from the surface with an intensity two times the bulk density.
In order to determine which properties are the driving force for the K ion to penetrate the inner Helmholtz layer, in contrast to the other cations considered, we performed a series of simulations with a negative charging of the electrode. In each simulation one cation was frozen either in the bulk solution or close to the surface; the surface absorption height was chosen according to the first adsorption peak in the density profile as reported in Table 1. Figure 2 reports the total number of coordinated water molecules within a given radius in bulk (a) and adsorbed on the surface (b). Within the first coordination sphere, which is 0.35 nm, the difference between the bulk and surface is minimal (< 0.1 water molecules) for all ions except K. The latter has on average 1.5 less coordinated water when it is adsorbed at the surface. These results are highlighted in Figures 2c and 2d, which show the coordination of a Na and K ion respectively at the charged surface. For Na, the closest 5 water molecules make up the first coordination shell, of which two clearly reside in between the ion and the surface. In contrast the K ion is directly adsorbed on the graphene surface with no intermediate water molecule. We note that even in the bulk K does not have a stable first solvation shell; unlike each of the other ions the K radial distribution function has a broad peak with an elongated tail which signifies that the final molecule in the coordination sphere is weakly bound at a moderately large distance from the ion. Consequently, there is a comparatively low energy requirement to remove such a molecule from the solvation shell with respect to the other ions investigated that permits the dehydration we observe at the surface. Behaviour which is consistent with results from investigations of K ion migration through biological channels. 37 Electrolyte

Figure 2: Computed radial distribution functions and integrated radial distribution functions of frozen cations in the bulk solution (a) and at the surface (b) for the various 1M concentration electrolytes considered in this work. Representative configurations of the Na (c) and K (d) ions frozen at the interface.
Direct exposure of the K ion to the surface leads to an induced polarization of the C atoms close to the ion resulting in strongly negative attractive Coulomb interactions. In Figure 3a, for a single example configuration, we plot the surface charges in two dimensions and mark the positions of K ions directly adsorbed to the graphene with crosses circled blue. As a guide to the eye, the surface charges are reported normalized to the homogenous fixed charge −3.4 336 ⁄ e (where −3.4 is the total charge of the graphene sheet made of 336 carbon atoms). In this example the surface Mulliken charges are significantly more negative in the surface region local to the ion.
To understand the extent to which this is a general property of the directly adsorbed ions, we also plot in Figure 3b all surface charges as a function of the distance from each adsorbed K ion across the 15 ns trajectory. As seen by the linearly increasing trend, up to 0.35 nm the polarization of the C atoms is very strongly influenced by the ion, this corresponds to nearest and immediately bonded atoms in the surface. Beyond 0.35 nm the computed Mulliken charges average out to the homogenous fixed charge value through interactions with other ions or water molecule dipoles. Our results demonstrate the importance of the local ion induced polarization in the stabilization of the direct adsorption configuration and the clear trend and distance dependence opens up for the possibility of an analytical description of the polarization energy to further expedite QM simulations. 38 We turn now to the electrochemical double layer integral capacitance, which is obtained from the system number density profile along the surface normal, ( ). First, the charge density, ( ), is calculated from the atom-resolved number densities multiplied by their partial charges from the classical force field. This includes all charged atoms (appropriately scaled 29 ), ions and electrodes. Poisson's equation relates the curvature of the electrostatic potential, Φ, to the charge density, which along one dimension reads, where 0 is the vacuum permittivity constant and is the specific electrolyte permittivity. Integration of Poisson's equation yields the electrostatic potential, Equation 2 is integrated from a position in the bulk electrolyte where the net charge density is zero, towards the quantum mechanically treated electrode. The bulk region has a constant electrostatic potential due to its net electroneutrality, which means that bulk potential can be defined as a reference for the calculation of the specific electrode capacitance, Φ Bulk = Φ ref . 39 The resultant potential drop across the interface ΔΦ = Φ electrode − Φ ref , which is depicted in Figure 4, at a given surface charge density, , can be used to obtain the electrochemical double layer integral capacitance, EDL , of the electrochemical double layer, EDL = ΔΦ . 3 The treatment of the electrolyte permittivity for this family of interfaces has been subject to recent scrutiny. 40 Standard protocol is to approximate using the dielectric constant of the electrolyte, yet it is well documented that truncation of the dipole field leads to a strong long-ranged dependence on the distance from the graphene surface. In fact, by fitting a polynomial function to the ( ) of SPC/E water from reference 41 Finney et al were able to extract the capacitance of graphite-NaCl solutions at various concentrations taking into account the changes in due to proximity to the electrode. 40 For our systems, given the different structuring of the water coordination shell shown in Figure 3b close to the electrode, the will be sensitive to the cation type. 9 Determination of ( ) requires prohibitively expensive s scale simulations of the interface, 41 thus we adopt an approach similar to Finney et al, fitting a function to the available data, scaled according to the dielectric constant we obtain for each electrolyte in the bulk phase.
Another facet of the electrochemical interface we can explore with the QMMD method, which is not accessible by purely classical approaches, is the contribution of the electrode electronic structure to the integral capacitance. Where conventional metallic electrodes, which have near infinite DOS at the Fermi-level, effectively screening the contribution of the electrode electronic structure, a graphene electrode is two-dimensional and semimetallic. The ensuing finite density of states (DOS) close to the Fermi level responds to an applied potential by populating the conduction (valence) band edge with negative electrons (positive holes). The rate of change of the voltage of the electrode is different from the potential difference, and is referred to as the quantum capacitance, Q , 42 which contributes to the total specific electrode capacitance as The differential Q can be directly obtained from the electrode DOS, which in turn can be obtained at the DFTB-level of theory at each iteration of the QMMD loop. Here is the energy relative to the Fermi level, ( ) is the DOS at a given energy and is the electron charge. Integration of the Q diff furnishes the integral quantum capacitance, 43 In practice, previous work has shown that the graphene DOS is resistant to changes in the presence of aqueous electrolytes, [44][45][46] as well as to charging of the sheet, 47,48 therefore a reasonable and commonly adopted approach is to use frozen bands; 49 the DFTB computed DOS and are plotted in Figure 4. For each of the systems considered here, the computed ΔΦ, EDL and , are reported in Table 1.  In our setup, simulations are carried out at constant surface charge, 50 with the distribution of the surface charge density modulated through feedback with the electrolyte. The computed EDL and therefore depend on the measured potential drop according to equations 3 and 4. First we consider the capacitance of the electrochemical double layers of the various ions, which are reported in Table 1. At constant negative surface charge, the EDL is roughly constant through the different monovalent ions, with LiCl having a slightly lower (~4 %) computed value. The EDL for the divalent ions are also found to be equal, however, these are significantly lower (~28 %) than the monovalent ions. On the other hand, at constant positive surface charge the EDL of all systems are found to be equivalent due to the independent adsorption of the Cl anion on the graphene in all cases (see Figure   1). Across all systems we find that the negatively biased electrode has a computed EDL nearly 5 times greater than the comparative positively biased electrode.
The reasons for the various similarities and differences between each of our simulations can be understood through examination of the ΔΦ. Factoring in the distant dependent ( ), we find that, for the same electrode surface charge density, the ions Li → K yield an ΔΦ ≈ −0.32 V. We note that LiCl solution has a slightly more negative ΔΦ, thence through the inverse relationship of Equation 3, the EDL is lower. What is most remarkable about this result is that, contrary to intuition, irrespective of the adsorption mechanism of the ion of the ion on the surface (K: direct adsorption through dehydration; Li, Na: indirect adsorption with complete 1 st hydration shell) the ensuing electrostatic potential drop arising across the double layer is the same over a 60 ns QMMD trajectory. In other words, the local structuring of the ions at the interface have a negligible effect on the ensuing double layer capacitance. Differences instead manifest in the simulations of the divalent cations; As shown in Figure 1, the concentration of the divalent ions in the electrochemical double layer is the same as that of the monovalent ions, without changes to the Cl concentration. As a result, Φ of the double layer is driven to more negative values and the ΔΦ are ≈ −0.1 V larger. This small change in ΔΦ mediated by the increased ionic charge leads to the significantly lower EDL values obtained.
We also checked the effect of doubling the ion concentration of the various electrolytes to 2M, apart from Li, we find that the computed electrostatic potential drop and integral capacitance do not change appreciably at higher concentration. Instead, Li is found to have a larger electrostatic potential drop due to the presence of more cations at the interface than anions. We speculate that the smaller Li cation can pack more ions at the surface and thereby decrease the capacitance.
Previous MD simulations have observed a different trend for solutions of NaCl, KCl, RbCl and CsCl, where the EDL decreased with ionic radius. 51 However, in reference 51 polarization effects are not accounted for and the dielectric permittivity of the aqueous solution was treated as 1. We note that we also observe a decreasing EDL by setting = 1. A very recent contribution from Dockal et al reports essentially the opposite trend, with EDL increasing with ionic radius; 52 in Dockal's work the polarizability of the surface is treated statically within the non-bonded Lennard-Jones parameters 19 and the is also set to 1. Yang et al addressed the dependence of EDL on in the Helmholtz layer by applying computed values of for each electrolyte. 9 In their work they computed a constant EDL for the cation series Li to Cs, Ca and for mixtures of the cations. The constant EDL is explained by the fact that the larger cations have both a larger dielectric constant and larger double layer thickness. According to the traditional Helmholtz theory, these two effects compete, resulting in a constant capacitance across the series, an effect we also observe in Figure 1. For our systems we obtained a constant double layer thickness and only marginally increasing dielectric constant for increasing ion sizes, moreover the divalent ions are found to have lower dielectric constant.
Turning to the electrode electronic structure and computed quantum capacitance, besides a rigid shift in the absolute position of the Fermi energy, the DOS is found to be consistent across each of the different electrolytes and surface charge density regimes. The fact that the energetic dispersion of the DOS does not change at different charge states indicates that the frozen bands approximation adopted for our simulations is valid. The computed DOS, reported in the top panel of Figure 5, shows that the DFTB approach captures the correct electronic features of graphene, namely the zero DOS and the linear dispersion of the band energies at the Fermi level. Integration of the DOS according to Equation 5 provides the differential quantum capacitance plotted in the bottom panel of Figure 5. The Q diff has a value of 2 µF cm -2 at zero potential and increases with positive and negative bias. These results are in agreement with other approaches to obtain the Q diff : for instance Xia et al employ a band theory model based on the two-dimensional electron gas (2DEG) and find that the minimum has a value of ~7 µF cm -2 . 53 Xia et al show that the 2DEG would have a zero Q diff without temperature induced broadening and impurity effects in the lattice. Interestingly, our results indicate only a minor asymmetry between the two branches of the differential capacitance, with the negative branch increasing slightly more rapidly, ( Figure 5, bottom panel). On this basis, it can be said that strong asymmetries in the experimentally measured total specific capacitance arise from the EDL component, i.e., preferential adsorption and structuring of either cations or anions on the surface. In other words, the Helmholtz capacitance for different ionic species cannot be considered constant.
Integration of Q diff over a potential window according to Equation 6, and represented pictorially in the bottom panel of Figure 5, yields the integral quantum capacitance. This can be combined with the EDL capacitance to calculate the total electrode specific capacitance as reported in Equation 4. As reported in Table 1, we find that the specific electrode capacitance (i.e. factoring the EDL and quantum effects) increases with ionic radius from Li to K and with ionic charge from +1 to +2 on the negatively charged electrode surface. On the positively charged surface the specific capacitance is fairly constant across ionic radii and charge. Although the absolute values of the specific integral and differential capacitance cannot be directly compared, our computed trends in can be said to be in agreement with the measurements of the basal plane capacitance of graphite electrodes. 10 Finally, we explore the dynamical properties of the ions at the interface. In comparison to the bulk liquid phase, water at the graphene-electrolyte interface was recently observed to have a higher selfdiffusion, 31,54 this is the result of a redistribution of the hydrogen bonding network at the interface increasing the local packing of water molecules on the surface. In addition to this, other dynamical properties such as the water residence time at the surface changes when explicit polarizability of the surface and molecules are factored into the model. 55 Motivated by these changes in the behaviour of liquid water at the graphene interface, here we explore the ensuing changes to the self-diffusion coefficient for the different cations in contact with the charged polarizable graphene electrode. This is achieved by means of the Einstein Mean Squared Displacement (MSD), 55 modified to measure the displacement of particles in the plane of the electrode surface. The cation diffusion coefficient, , is calculated from the slope of the MSD over time, where ( ) is the particle survival probability, this is nothing more than the probability that a given cation will remain within a slab of thickness Δ in the time interval .
First we consider the self-diffusion of the water and ions in the bulk phase, these have previously been benchmarked for our simulation setup elsewhere 29,57 and help to validate our investigation on the interface. Far away from the electrode we recover the bulk diffusion coefficients of Li (0.083), Na (0.103), K (0.166), Mg and Ca (0.058) to within 0.01 nm 2 ps -1 . These results match trends observed for ion conductivities in the bulk: namely that conductivity, thence ion diffusion, increases from Li to K and decreases moving from an ion with oxidation state +1 to +2. 58 For analysis of the interfacial diffusion, we partition the simulation box in to 0.75 nm slabs, this encompasses the local structuring of each ion at the surface. We find that the rate of diffusion increases for each of the ions as reported in Table 2. This is consistent with the picture for water (in contact with a single graphene sheet), 31,54 where the reorganization of the hydrogen bonding network into pentagonal rings at the surface results in a higher local packing of the water molecules compared to the bulk. While polarization effects are known to increase the residence time for water molecules on the surface, the diffusion coefficient is only minimally affected. 55 Therefore, given the large disparity in the number of ions vs water molecules at the surface, its most likely that the enhanced water molecule diffusion at the surface induces a similar effect on the ions in the outer Helmholtz layer. When the electrode is charged, we see very little change to the self-diffusion of cations at the surface. Analysis of the H-bonding network of the water molecules seems show that the H-bonding network is only negligibly affected by the charge on the electrode. 52 Consequently, the diffusion of the water thence the diffusion of the extremely low number of ions in that water is also only negligibly affected.  The changes to moving from bulk towards the surface could also be the result of changes to the local concentrations of ions within EDL. The typical trend is that decreases with concentration, 29 while in our simulations the concentration of ions close to the surface (7.5 Å) is larger than in the bulk, as depicted in Figure 1. However, we can rule out the contribution of concentration since upon charging, the concentration of cations at the interface increases by at least a factor of 2, with no obvious change to their diffusion. That is, except for the K ion, which is significantly slowed with respect to the uncharged case. In the charged system, the concentration of K ions in the 1 st layer (7.5 Å) is close to three times the uncharged case, and the number of ions that penetrate the inner Helmholtz layer is also greatly increased, Figure 1. However, the change to K is also likely a surface effect: as we showed in Figure 3a and b, K ions in the inner Helmholtz layer strongly polarize the charged surface, this strong polarization-enhanced Coulomb interaction binds the ions to the surface atoms resulting in decreased in-plane mobility. The slowing of specific ions, based on their ability to bind to the surface, thence their hydration free energies, raise the interesting possibility for charged graphene electrodes, nanochannels and nanoflakes to act as capture devices for specific ions in the separation of mixed ionic solutions. We note, this dynamical effect cannot be observed in classical models that do not take explicit polarization into account, nor those which restrict fluctuations of the surface charge density to the locality of the atom since the requirement for the accumulation of charge at the adsorption site cannot be met. Moreover, the difference between the behaviour of K and the remaining ions Li, Na, Mg and Ca can only emerge in the case where the dynamics of water and ions are be treated on an equal footing. In this respect, where fully ab-initio methods can prove too costly, the lack of a solvation shell in implicit solvent models contributes to the failure towards a failure to distinguish between ion adsorption in the inner and outer Helmholtz plane.
We can also compare the dynamical behaviour of the ions at the interface to recent experiments on the conductivity of ions in nanochannels, 58 i.e those with an interlayer separation equivalent in size to the hydrated ion. In those experiments, the mobility of Li and Na was found to increase moving from a bulk electrolyte to channel, K was found to have a roughly equivalent mobility while Mg and Ca both slowed. The change in behaviour is linked to the diameter of the hydrated ion, its hydration free energy and the ability for its hydration shell to distort or deplete. As we shown in Figure 3a and b, we find that the 1 st solvation shell of the ion is unaffected by the presence of the surface except in the case of K. We also checked the 2 nd solvation shell, which further confirms that only the solvation shell of the K is notably affected via adsorption on the surface. Therefore, taking into account our results and those of reference 58 , we tentatively suggest that the increased mobility of the Li and Na observed in experiment is the result of exposure of the ion and its solvation shell to the graphene surface. The equivalence between the bulk and confined mobility of K results from a cancellation of the speed up due to its exposure to the surface and the direct adsorption of the ion polarizing either (or both) of the channel surfaces. Finally the reduction in mobility of the Mg and Ca, which have both the largest ionic charges and solvated diameters, due to the forced dehydration of the ion entering the nanochannel and the highly charged ion strongly polarizing the surfaces and slowing down the diffusion dynamics. Clearly, whilst these insights provide promising explanations of the observed trends, further investigation, and validation of the behaviour of the ions within subnanometer cavities is required.
In summary, our QMMD simulations provide a new take on the contribution of specific ion effects to the thermodynamic and kinetic properties of aqueous graphene-electrolyte supercapacitors. We find that the coupling and interplay between dehydration and induced surface polarization dictates whether ions can penetrate the inner Helmholtz plane and directly adsorb on the surface. However, changes to the electrode specific capacitance with ionic radii are only obtained when the quantum capacitance is included. Regardless, complete discharge of the surface can be obtained by oppositely charging the electrode. Our calculations revealed that any asymmetries between the positive and negative branches of the measured capacitance originate in the EDL capacitance. This has implications for the modelling of a constant Helmholtz capacitance between solvated anions and cations. Finally, we found that the diffusion of all ions at the interface is greater than their bulk values, in line with results for water. K is the only ion in our study to slow down upon charging of the electrode, this is due to the induced polarization of the graphene surface and highlights the important role surface effects have on ion mobilities on graphene, graphite and within graphitic nanochannels. Our results can inform on the future design of electrode-electrolyte coupling and provide a basis for analytical treatment of the quantum mechanical polarization.