Microscopic Insight to Nonlinear Voltage Dependence of Charge in Carbon-Ionic Liquid Supercapacitors

The impact of cell voltage on the capacitance of practical electrochemical supercapacitors is a phenomenon observed experimentally, which lacks a solid theoretical explanation. Herein, we provide combined experimental and molecular dynamics investigation of the relation between voltage and capacitance. We have studied this relation in supercapacitor cells comprising of activated carbon material as the active electrode material, and neat ionic liquids (ILs), and a mixture of ILs as the electrolyte. It has been observed that the increase of accumulative charge impacts the conformation and packing of the cations in the anode, which determines its nonlinear behavior with increasing voltage. It has also been shown that for the mixture IL with two types of cations, the contribution of each type of cation to the overall capacitance is highly dependent on the different pore sizes in the system. The smaller tetramethylammonium (TMA) shows tendency for more efficient adsorption in the mesopores, while 1-Ethyl-3-methylimidazolium (EMIM) is found to be present almost exclusively in the micropores where TMA is present in small quantities. Such microscopic insights from computer simulations of the molecular phenomena affecting the overall performance in supercapacitors can help to design more efficient electrolytes and devices.


Introduction
Development of new and more effective electrochemical energy storage, such as high energy and power electrochemical supercapacitors, is in high demand with the increase of renewable energy produced globally [1][2][3]. Ionic liquids are promising electrolytes of choice for the electrochemical energy storage application due to their large operational voltage. Improving the capacitance of high-voltage supercapacitor in IL electrolyte is an efficient way to achieve higher energy and power density. In different studies, the increase of the applied voltage leads to an increase of the capacitance in general, although the experimental setups differ from each other [4,5]. The effect of the applied voltage on the increasing capacitance can be up to 50% for activated carbon (AC) electrodes and up to four times for singlewall carbon nanotube (SWCNT) [4,5]. Both linear and nonlinear dependence of capacitance with the voltage has been observed in these studies.
For an ideal electric double layer capacitor [4], the charge increases linearly with the increase of the applied voltage, and therefore, the capacitance should remain constant. However, many parameters, such as the complicated porous structure [5] and defective carbon surface, could lead to a nonlinear relationship between voltage and charge. Recently, we have shown that the more effective ionic packing in the confined pores is a key parameter for increasing the capacitance [6]. Moreover, the interactions between different cations used as a mixture in the same supercapacitor cell, such as 1-Ethyl-3-methylimidazolium (EMIM + ) and tetramethylammonium (TMA + ), can also alter the packing density in the porous electrodes and, in effect, alter the capacitance [7]. The orientation of nonspherically shaped ions (EMIM + ) changes with the applied voltage, which can alter the packing density of the ion on the electrode surface and the capacitance per surface area [8]. Hence, fundamental insight and understanding of the detailed interactions between the electrode surface and the adsorbed ions at different applied potentials are essential in order to rationally tailor the porous electrode material and IL composition for improved SC prototypes.
Molecular simulations are often used to provide qualitative and quantitative insights related to molecular scale phenomena to explain the behavior of complex physical and chemical systems. In particular, molecular dynamics (MD) simulations [9][10][11][12][13][14][15][16] have been successfully applied to study and understand the molecular phenomena occurring in supercapacitor systems. The electrodes (both anode and cathode) in SC have been modelled by a single sheet of graphene [9][10][11][12], carbon nanotube [13][14][15], or complex nanoporous carbon material based on experimental observations [16]. The distribution of the ions and the properties of the electrical double layer on the anode and cathode sides have been a common point of discussion in the simulation studies. In all of these systems, neat ionic liquids (IL) [9][10][11][12][13][14][15][16] or IL in a mixture with organic solvent (acetonitrile) [11] have been used as the electrolyte. Shim et al. [11] have also evaluated the effect of acetonitrile mixed with IL on the performance of the supercapacitor, comprising of graphene sheets, and have shown that the capacitance of the pure IL (EMIM + BF − 4 ) is higher than the capacitance of a mixture with acetonitrile. However, the interactions between different cations in the electrolyte mixture in the same supercapacitor cell have not been investigated thoroughly [7].
The structure and density of the ions near the electrodes (both positively and negatively charged) represented as a single graphene layer have been extensively investigated [9][10][11]. The high capacitance of the IL supercapacitors is correlated with the strong and effective screening of the charged electrode surface by the IL [12]. Nanoporous supercapacitors have been modeled by a carbon nanotube (CNT) [13][14][15] or carbon material having complex geometry based on experimental pore size distributions [12]. Shim and Kim [13] have studied the behavior of ions inside the pores and the capacitance for various diameters of the CNT for the EMIM + BF − 4 ionic liquid. Their results [13] suggested that the capacitance of a CNT supercapacitor with a diameter in the range of 0.9-2 nm increases with decreasing CNT diameter. That observation was associated with the different internal solvation exhibited by ions inside the pores. It was also shown that the lack of internal solvation for pore diameters below 0.8 nm significantly decreases the capacitance. Although many MD studies [9][10][11][12][13][14][15][16] provide valuable information for understanding the molecular phenomena involved in SC systems, they have not focused on the relations between capacitance, charge, and voltage in a SC.
In the current study, both atomistic MD simulations and experiments have been performed to reveal the voltagedependent charge density relation in porous carbon-ionic liquid supercapacitors. The MD simulations were performed in order to obtain the densities and spatial distributions of IL ions inside a negatively charged porous electrode, modelled by a single wall carbon nanotube. The densities and spatial orientation of the ions inside the pores were calculated for two different pores sizes, four different surface charge distributions, and three different IL compositions. We have determined how the spatial orientations and the densities inside the pores influence the relation between volumetric charge density and voltage in various pores and mixtures of ions in electrolytes. The simulations are shown to be qualitatively in agreement with the experimental data. The simulation setup presented here, although exploring an idealistic case study of a real supercapacitor porous electrode, is suitable for studying the ion orientations and densities and their impact on the capacitance in an ionic mixture electrolyte.

Methods and Models
2.1. Molecular Dynamics Simulations. A simulation box of 500 pairs of 1-Ethyl-3-methylimidazolium (EMIM + ) and tetrafluoroborate (BF − 4 ) was initially created. In addition, 110 and 205 pairs of tetramethylammonium (TMA + ) and BF − 4 were added to mimic the IL mixture composition in the experiments. In total, three simulation boxes were created with 0 wt% TMA + (System 1), 15 wt% TMA + (System 2), and 25 wt% TMA + (System 3) diluted in the EMIM + BF − 4 ionic liquid (IL). The simulation boxes were initially equilibrated in the isothermal-isobaric statistical (NPT) ensemble at 293.15 K and 1 bar of pressure, and the densities of the IL were calculated. The density of the neat IL (0 wt% TMA + ) was calculated to be 1293 kg/m 3 which is in very good agreement with the experimental value of 1297 kg/m 3 [17].
Subsequently, we have used these bulk simulation boxes to create the porous electrode simulations. The porous electrode was created by joining a carbon nanotube with two graphene layers as depicted in Figures 1(c)-1(f). Two of the bulk simulation boxes were placed on both sides of the graphene layers. The porous electrode and simulation snapshots were created with the Visual Molecular Dynamics (VMD) package [18]. In this respect, an idealistic case of a true pore in the real porous electrode has been built, which will enable us to examine the spatial and density distributions of the ILs inside micro-and mesopores in order to construct an abstract model of cation adsorption for the SC anode system described in Reference [6]. Two different carbon nanotubes with diameters of 1.085 nm and 2.169 nm have been used to represent a general micro-and mesopore, respectively. The length of the pore was 4.919 nm and was chosen as a compromise between computational efficiency and sufficient space for the ions to migrate inside. The two graphene layers with dimensions of 4:9 × 4:9 nm were placed on both ends of the nanotube. The pore was then left to be filled with the ions from the bulk, which was the first equilibration for 500 ps in the NPT ensemble. A second equilibration also in the NPT ensemble was then carried out for 3 ns. This was done in order to achieve the correct bulk value of the density outside the pore when some of the ions migrate into the pore. The pressure coupling, maintained by the Berendsen barostat [19] at 1 bar, was applied semi-isotropically, where the xy-dimensions of the simulation box were kept constant and only the z-dimension was scaled. The density was then calculated again after the 3 ns equilibration to assure that the correct bulk density was achieved. After that, a 5 ns equilibration in the canonical statistical ensemble (NVT) was carried out, which was followed by a 10 ns production run also in the NVT ensemble. The coordinates 2 Energy Material Advances of the molecules during the production run were saved every 5 ps. The nanotubes were negatively charged in order to mimic the applied voltage on the anode in the experimental setup. Four different surface charge distributions (-0.23 e/nm 2 , -0.53 e/nm 2 , -1.00 e/nm 2 , and -2.00 e/nm 2 ) were modeled by assigning negative point charges to the nanotube atoms. These surface charge distributions should approximately correspond to -1.0 V, -2.0 V, -4.5 V, and -9.0 V applied to the anode, respectively [12]. The numbers of the cations and anions for the different surface charge concentrations of the pore were adjusted in order to keep the overall electroneutrality of the systems. The total charge on the micropore was -4e − , -9e − , -18e − , and -33e − corresponding to the surface charge distributions of -0.23 e/nm 2 , -0.53 e/nm 2 , -1.00 e/nm 2 , and -2.00 e/nm 2 , respectively. Therefore, we removed the respective number of BF − 4 ions to keep the systems electroneutral. In the mesopore systems, we removed 8, 18, 36, and 67 anions, respectively. It should be noted that this is a standard procedure found in many simulations [11,13,15].
The EMIM + , BF − 4 , and TMA + ions were modelled in full atomistic resolution by using the force field parameters given by Tsuzuki et al. [20]. These force field parameters are slightly modified OPLS-AA parameters for ionic liquids based on those developed by Lopes et al. [21,22] with the addition of parameters for the TMA + cation. The geometrical combination rule has been used for the nonbonded interactions which are in accordance with the OPLS-AA force field. The temperature was set to 293.15 K and maintained by the velocity rescaling thermostat [23] with a time constant of 100 fs. The time step for integration of the equations of motion was set to 1 fs. The electrostatic interactions have been modelled by the Particle-Mesh Ewald summations with a cut-off distance of 1.2 nm. The MD simulations were performed with the GROMACS 5.0.4 simulation package [24], and periodic boundary conditions were applied to all three spatial directions.
The densities, spatial distributions, projected surface area per molecule, and orientation angles of the ions inside the pores were calculated. A molecule is counted to be inside the pore when its center of mass (COM) is inside the nanotube, irrespective of the time spent inside the pore. The projected surface area per molecule was calculated for the TMA + and EMIM + cations, using the Van der Waals molecular surface area projected on a flat surface. Spatial distributions in the radial direction of the nanopores were calculated as follows: the COM of every molecule was calculated first, and then their positions were sorted in shells inside the pores, starting from the pore center, where the shell thickness was 0.1 Å. In that way, the probability of finding a 3 Energy Material Advances given molecular type in a given shell in the radial direction was obtained, as well as the distance to the nanopore wall. The orientation angles only for the EMIM + cations were calculated, since the other ions are symmetric. Moreover, EMIM + is more flexible compared to TMA + and can adopt different conformations depending on the environment. Therefore, its orientation inside the pores can be influenced by the surface charge.

Experimental
Methods. Three microporous carbon and three mesopore-rich activated carbon electrodes were studied in this work, and all materials' preparation method has been described in detail in our previously published work [7]. Three mesopore-rich carbons were prepared by chemically activating carbonized polyaniline with different carbonization and activation temperatures. Three microporous carbon are commercial activated carbon YP-50, carbide-derived carbon, and activated carbon nanospheres, respectively. Detailed pore characteristics and surface area information of all carbon materials are listed in Table 1. For more detailed information of the material preparation, please see Reference [7].
The carbon materials were separately milled with 5 wt% Polytetrafluoroethylene (PTFE) solution (Sigma-Aldrich, containing 8.4 wt% PTFE) in a small amount of ethanol. Afterwards, the mixture of carbon and PTFE was kneaded and rolled into a film and then cut into a disk-shaped electrode with a diameter of 8 mm and mass of 2.5 to 4 mg. The disk-shaped dots were dried in vacuumed oven overnight before being used as active electrode materials for the supercapacitors.
Symmetric configuration of the SCs was used in this study, with two carbon electrodes with the same mass as positive and negative electrodes, respectively. A carboncoated aluminum foil was used as a current collector, while a 25 μm thick polypropylene membrane (Celgard 3501) was used as a separator. All cells were fabricated in a glovebox (Vacuum Technologies Inc.) with O 2 < 0:1 ppm and H 2 O < 0:1 ppm. After adding 1-Ethyl-3-methylimidazolium tetrafluoroborate (EMIM + BF − 4 ) or 15 and 25 wt% of tetramethylammonium tetrafluoroborate (TMA + BF − 4 ) in EMIM + BF − 4 as electrolyte onto the electrode materials, the cells were sealed and stabilized overnight before testing. The capacitance of a single electrode was calculated from the discharge curve of Galvanostatic charge/discharge based on where I is the current used during discharging, m total is the overall mass of the active materials on both electrodes, Δt is the discharge time, and V cell is the voltage change during discharging (excluding the voltage drop).

Charge/Voltage Relation in Micropores.
The results obtained from the MD simulations and the experiments for the voltage/charge relation in the micropores are presented in this section, while the voltage/charge relation in the mesopores is presented in the following section. The capacitance is normalized by the micropore volume of each material to calculate the contribution of capacitance per pore volume-C vol (F/cm 3 ). Further, the volumetric charge-ρ (C/cm 3 )-is estimated by the following formula: where V is the voltage. Similarly, ρ from the MD simulations has been calculated in C/cm 3 , where the average number of cations inside the pore per unit of volume has been calculated from the production run. Furthermore, ρ has been normalized with respect to the lowest TMA + concentration case and lowest voltage (surface charge density in the simulations) (ρ 0 ), for both experimental/simulation results, and referred to as normalized charge density (ρ/ρ 0 ). In that way, we can compare the trends obtained from the simulations with the experimental results. Figure 1 presents the voltage/normalized charge relation in the micropores obtained from the experiments (Figure 1(a)) and simulations (Figure 1(b)), as well as simulation snapshots for the different surface charge densities employed in the simulations. The normalized volumetric charge increases when the voltage (or surface charge density) increases, correlating to the filling of the pore with ions. The slope of the curves gradually decreases, which can be observed from both the experimental and simulation result, indicating the capacitance gradually decreases as the potential increases. Figure 1(a) shows that the increase of TMA + concentration in the experimental system decreases the capacitance in the micropore. On the other hand, the simulation results in Figure 1(b) demonstrate a negligible increase of the normalized charge with the increase of the TMA + concentration, and it is explained in the next paragraph. Table 2 summarizes the results for the average number of ions in the pores ðN E , N B , N T Þ and their respective number    Figures 1(c)-1(f) and Table 2 show that mainly EMIM + cations are present inside the micropore when the surface charge density is low. As the surface charge density increases beyond -1.00 e/nm 2 , TMA + starts to present in minimal quantities in the micropore. The additional TMA + molecules, for surface charge density of -1.00 e/nm 2 and -2.00 e/nm 2 , contribute to the increase of the charge density since the surface area (0.4 nm 2 /molecule, Table 2) occupied by them is smaller, compared to EMIM + (0.56 nm 2 /molecule, Table 2), and therefore, on average, more molecules can enter the pores. For surface charge density of -0.50 e/nm 2 , it can be seen from Figure 1(b) that the volumetric charge density is the same for all three systems studied. This can be explained with the fact that for all systems with this surface charge (-0.50 e/nm 2 ), the same number (7) of EMIM + molecules is found inside the micropores, irrespective of the TMA + concentration. For the lowest surface charge density, also some BF − 4 molecules are found in the micropore, and they are slightly more on average for systems 2 and 3 (≈1.2) compared with system 1 (≈0.2). The effect of this is that more EMIM + molecules are found in system 2 and 3 (≈5) compared with system 1 (≈4), as seen from Table 2. Thus, the volumetric charge density in systems 2 and 3 is higher compared with that in system 1. Therefore, the increase in the surface charge density for systems 2 and 3 (except for -0.50 e/nm 2 , which is identical for all systems) in the simulations is explained with the increase of BF − 4 and TMA + molecules present in the micropores in these two systems. For the two lowest surface charge densities, no TMA + was found in the micropore, consistent with our previous study [7]. Meanwhile, the TMA + with a more concentrated charge should have a stronger attraction with the electrode (evident also from the radial distributions presented in Section 3.3, where TMA + is found closer and more tightly packed than EMIM + to the pore wall); hence, the TMA + ion transport is more resistive compared with EMIM + . With such a confined space, the presence of the TMA + can easily block the pore, leading to a lower capacitance. This explains the noticeable capacitance decreases when TMA + is added based on the experimental result. Especially when the voltage increase, the excess of the TMA + ion leads to more capacitance loss as shown in Figure 1(a).

Charge/Voltage Relation in Mesopores.
The relation between normalized charge and voltage in the mesopores is presented in Figure 2. Both experiments (Figure 2(a)) and simulations (Figure 2(b)) clearly show that the capacitance increases with the addition of TMA + . The slope of the curve becomes smaller experimentally, indicating capacitance decreases with the increase of the applied voltage. The smaller slope at the high voltage agrees with the observation in the MD simulation for the 0 wt% and 15 wt% of the TMA + . Larger ion transport resistance (pore blockage due to stronger interfacial forces, as explained in Section 3.1) at larger potential, may be one explanation of decreasing capacitance. However, an abnormal capacitance increase in the simulations (Figure 2(b)) is observed at large applied voltage for 25 wt% of the TMA + . This is probably due to an abrupt increase of the ratio of the TMA + (n T /ðn E + × n B + × n T + Þ) from 15% to 22.4% when the surface charge density increases from -1.00 e/nm 2 to -2.00 e/nm 2 .
From the comparison of Figure 1(b) with Figure 2(b), it can be noted that the increase of the normalized charge with voltage is more than twice as effective in the micropores than in the mesopores based on the simulation result. This indicates that the capacitance per pore volume is much larger in the micropore than in the mesopore, or in other words, the utilization of the pore volume by the ions is more efficient in the micropores. This agrees with the experimental results and the literature that the normalized capacitance is higher in the micropores [6,25].
The surface charge effects on the flexibility of EMIM + in the mesopore are expected to be less pronounced due to the screening of the charges by the presence of the BF − 4 ions in the pore. The presence of BF − 4 in the mesopore is significant in the lower surface charge cases, since the negatively charged surface is screened by the cations. However, with the increase of the surface charge, the amount of BF − 4 progressively lowers, as evident from Figures 2(c)-2(f) and Table 2, which leaves more space for cation adsorption. The effects of the surface charge on the density, spatial orientations, and distributions of the ions in both pores are examined in more detail in the next section.
To further study how the voltage affects the capacitance, we calculated the differential capacitance (C d ) as a function of voltage from the MD simulations. The differential capacitance was calculated similar to Reference [26], and details of the calculations can be found in the Supporting Information. It is evident from Figure 3 that with the increase of the voltage the capacitance decreases, as it was also explained above. This is also consistent with previous MD simulations of similar systems [15]. With increasing TMA + concentration, the capacitance increases for both micro-and mesopores in the simulations, which is also apparent from the voltage/normalized charge relation presented in Figures 1(b) and 2(b). Thus, it can be concluded that the addition of TMA + should, in the ideal case scenario explored by the MD simulation, increase the capacitance irrespective of the pore size. However, as shown in Figure 1(a), the addition of TMA + decreases the charge density inside the micropore leading to the decrease of the capacitance, which is attributed to more resistive TMA + ion transport compared to EMIM + .

Voltage Effects on the Ion Density, Spatial Distribution, and Orientations in the Micro-and Mesopores.
More information can be obtained by examining the number density of the different ions and particularly the number density of the cations in the pores. At the lowest surface charge density, the number density of the cations in the micropore is significantly lower than that in the mesopore. However, with the increase of the surface charge, the difference between the densities of the cations inside the two pores becomes insignificant for all the systems. Therefore, the packing of the molecules inside the micropore becomes more effective with the increase of the surface charge. The higher cation density in the pores for the elevated voltage would imply higher charge density. It should be also noted that the increase of the charge 6 Energy Material Advances density in the micropores comes from the increase of the EMIM + , since TMA + barely presents in the micropore. However, the opposite is true for the mesopores where the increase of the charge density in the mesopores comes from the increased presence of TMA + which also possesses smaller surface area than EMIM + . Therefore, a higher density of the  7 Energy Material Advances cations can also be achieved with the increase of TMA + concentration in mesopore systems.
That result is also in agreement with the results shown in Figures 2(a) and 3(b), where the mesopore capacitance increases with the TMA + concentration. It is also apparent from Figure 1(a) that the capacitance in the micropores even decreases with the increase of TMA + concentration. That moreover supports the idea of the contribution of TMA + to the increase of capacitance only in the mesopores and of EMIM + in the micropores. The increase of capacitance in these systems was described also in our previous study [7] and attributed to the selective behavior of the cations towards the micro-and mesopores. The results presented here strongly support our previous observations concerning the selectivity. Additionally, the results presented here suggest that EMIM + and TMA + show preferential adsorption in the micropore and mesopore, respectively. The selectivity behavior of the cations is further discussed in the next section. Figures 4 and 5 present the positions of the COM of the different molecules inside the micro-and mesopores (system 3) along the radial direction of the pores, respectively; here, zero on the x-axis represents the middle of the pore, and the probabilities are calculated in the radial direction with respect to the centerline of the pore. The dashed line in the figures represents the pore wall. When the surface charge is low (-0.23 e/nm 2 ) both the cations and anion can be found near the pore wall. The surface charge density is not high enough to repel the anion and could be screened by the cations. When the surface charge density increases, both the EMIM + and TMA + are brought closer to the pore wall, and TMA + is closer to the pore wall than EMIM + for both pore sizes, as seen from Figures 4(c), 4(d), 5(c), and 5(d). Inside the mesopore (Figure 5), we can further observe the formation of three distinctive layers and these layers become more prominent with the increase of the surface charge of the pore. Figure 6 presents the orientation angle distributions of EMIM + inside the micropore system 1 (Figure 6(a)), micropore system 3 ( Figure 6(c)), mesopore system 1 (Figure 6(c)), and mesopore system 3 ( Figure 6(d)). The orientation angle is calculated with respect to the axial direction (z-dimension) of the nanotube, where 0°is parallel to the z -dimension and 90°angle is perpendicular to the z-dimension, as depicted in Figure 6(e), where also the definition of the angle is presented. The orientation angle is defined as the angle between the vector connecting the two carbon atoms located next to the imidazolium ring and the axial direction of the nanotube. The surface charge of the nanotube affects the orientation of EMIM + significantly. When the surface charge is low, below 1.00 e/nm 2 , the orientation has a broader peak centered between 40°and 60°, whereas, 8 Energy Material Advances for surface charge of 1.00 e/nm 2 and above, several peaks appear where the most dominant peak is close to 90°. The peaks are also narrower compared to the peaks with the lower surface charge. This indicates also that EMIM + is almost locked in this configuration and a rearrangement of the molecules inside the micropore is difficult. Thus, the high surface charge has a solidifying effect on the molecules; i.e., they are close to a solid-state form as the surface charge increases. As seen also from Figure 6(g), EMIM + is quite flexible, allowing it to obtain configurations which follow the curvature of the micropore. Considering that our simulations do not include polarization effects, the dependence of EMIM + configuration of the surface charge is clearly seen. The orientation of EMIM + in the mesopores (Figures 6(c) and 6(d)) is slightly different, because of the presence of BF − 4 and higher concentrations of TMA + . It seems that the presence of TMA + shifts the orientation angle of EMIM + towards 90°. Also, with the increase of the surface charge, the angle shifts to 90°in the mesopore, but to a lesser extent.

Selectivity
Behavior and Free Energy of Adsorption. In our previous study, we have observed a selective behavior of the two cations towards entering the micro and mesopores [7]. It was shown that TMA + was able to enter only the mesopores, except when a higher surface charge was applied on the micropore. That was attributed to the stronger interactions of TMA + with BF − 4 in the bulk, which promotes more difficult dissociation of the cation-anion complex. In this section, we further explore the selectivity behavior by calculating the radial distribution functions (RDF) of EMIM + -BF − 4 and TMA + -BF − 4 (Figure 7) in the bulk and the respective free energy of adsorption into the pores of EMIM + and TMA + (Figure 8).
The RDFs presented in Figure 7(a) show that the interaction between TMA + and BF − 4 is much stronger than the interaction between EMIM + and BF − 4 . It is also apparent, from the decreasing of the first peak in the RDF for the different systems, that the interaction between EMIM + and BF − 4 becomes slightly weaker with the addition of TMA + . The first peak of 10 Energy Material Advances the RDF also reveals that the distance between TMA + and BF − 4 is around 0.4 nm, while the distance between EMIM + and BF − 4 is around 0.5 nm. This peak is also narrower for TMA + -BF − 4 compared with EMIM + -BF − 4 , which indicates the formation of more compact and more tightly bound first coordination shell around TMA + . As seen from Figure 7(b), TMA + is able to coordinate around 5 BF − 4 molecules while EMIM + coordinates around 7 BF − 4 molecules, which also indicates more loosely bound first coordination shell in the case of EMIM + . Therefore, the stronger interaction between TMA + and BF − 4 is likely the dominant factor that dictates the selective adsorption of the cations. TMA + can not enter the micropore, because of the strong interaction with BF − 4 , and also, their complex is large enough to prevent the adsorption in the micropore. However, when the surface charge increases beyond some threshold, the Coulombic interactions of TMA + with the nanotube become stronger and are able to break the complex with BF − 4 , and some TMA + cations are able to adsorb into the micropore. The adsorption in the mesopore, however, has a different mechanism. The dimensions of the mesopore are large enough to accommodate the complex of TMA + and BF − 4 which most likely are dragged together into the pore.
Here, we have also calculated the energy required for a cation to enter the micropore, which effectively is the energy penalty to move from the bulk to the pore. To wit, we have calculated the free energy needed for a cation to enter the micropore which is already filled by other ions from the electrolyte, thus obtaining the insertion free energy of a cation in the pore. This has been done by umbrella sampling simulations [27][28][29][30], and the free energy has been calculated by the Weighted Histogram Analysis Method (WHAM) [31,32]. The umbrella sampling procedure has been applied for system 3, for surface charge of -0.53 e/nm 2 and -2.00 e/nm 2 and for both TMA + and EMIM + cations. Either TMA + or EMIM + cations have been pulled from the bulk of the box in the direction towards the pore (z-dimension in the simu-lation box). The reaction coordinate for calculating the insertion free energy has been defined from the initial position of the chosen ion to the center of the pore. On average, 40 umbrella windows separated by 0.1 nm have been used to calculate the insertion free energy.
We have calculated the following free energies (Table 3): EMIM + (-0.53 e/nm 2 -25 kcal/mol), EMIM + (-2.00 e/nm 2 -45 kcal/mol), TMA + (-0.53 e/nm 2 -60 kcal/mol), and TMA + (-2.00 e/nm 2 -50 kcal/mol). The respective free energy profiles have been presented in Figures 8(a) and 8(b) for surface charge density of -0.53 e/nm 2 and -2.00 e/nm 2 , respectively. It is important to note, in the simulations with a surface charge density of -0.53 e/nm 2 , no TMA + molecules are observed to enter the micropore, and for -2.00 e/nm 2 , some TMA + molecules have entered the micropore during the nonumbrella sampling simulations presented in Section 3.1. Therefore, we have chosen these two cases as examples when TMA + is able or not able to enter the micropore and have calculated the energy needed for the molecules to enter the micropore. That should provide information about which molecule is more energetically favorable to enter the micropore. At the lower surface charge density of -0.53 e/nm 2 in Figure 8(a), the energy of TMA + near the beginning of the micropore is higher than the energy of EMIM + indicating that it is favorable for EMIM + to enter the micropore. At the higher surface charge density (-2.00 e/nm 2 ), the energy of the two molecules is nearly the same at the place where the micropore begins. The insertion of EMIM + in this case is still more favorable, but the energy difference between EMIM + and TMA + is only 5 kcal/mol. As it was already shown in Table 2 and Figures 1  and 4, some TMA + molecules are able to enter the micropore for the higher surface charge densities. It should be noted that the free energy of adsorption of TMA + decreases with the increase of the surface charge density while the opposite is observed for EMIM + . The former could be explained with the stronger attraction of TMA + to the negatively charged pore wall than EMIM + because of the higher charge density

11
Energy Material Advances of the smaller ion (TMA + ). Smaller ions with higher charge density interact stronger than large ions. Stronger attraction of TMA + to the pore wall is also evident from the radial distribution presented in Figures 4 and 5, where TMA + is located closer to the wall. Therefore, the increase of the surface charge density favors the adsorption of the smaller cation. The latter, the increase of EMIM + adsorption free energy, we speculate arises from competition of EMIM + for adsorption with TMA + .
From the umbrella sampling simulations, it can be also seen that the adsorption of EMIM + and TMA + is an endothermic process. We speculate that this is due to the formation of an electric double layer near the graphene surface as can be seen from the density plot presented in Figure 8(c) for the case of -2.00 e/nm 2 surface charge. Therefore, the cation needs energy to penetrate a region with an enhanced density, which makes the process endothermic. Figures 8(d) and 8(e) also present the trajectory of EMIM + and TMA + , respectively, to migrate from the bulk into the pore. From these trajectories, it can be observed that TMA + spends relatively more time near the graphene surface than EMIM + . EMIM + enters the micropore directly probably because it can also more easily change its conformation in the denser double  layer near the graphene wall outside the pore, than the more rigid TMA + . Another possible factor that can hinder the adsorption of TMA + is its stronger complexation with BF − 4 molecules (seen from the RDFs presented in Figure 7(a)) in the double layer region near the pore.

Conclusions
The relation between the applied voltage and volumetric charge density in a supercapacitor system has been investigated. The experimental results showing a decrease in capacitance with an increase of the applied voltage are compared with molecular simulations of a model electrode at corresponding surface charge densities. The simulation results are able to predict the experimental trends of the nonlinear relationship between the volumetric charge density and the applied voltage and the decrease of capacitance with increase of voltage. The molecular simulations give qualitative and quantitative insight of the ion packing inside the pores and can help with the understanding of the nature of the change of capacitance and charge with the increase of the applied voltage observed experimentally. Namely, it is shown that the different pore sizes contribute differently to the capacitance and total charge adsorbed in the pores and demonstrate the selective absorption of cations in the pores. EMIM + is mainly occupying the micropores, where only at high voltages TMA + is able to enter them, which is not apparent when only the size of the ions is taken into account. The more flexible EMIM + cation can easily enter the smaller pores and adapt better to the pore geometry and curvature, which leads to a more effective packing of the smaller volume of these pores. This is also suggested from the umbrella sampling simulations where lower energy is needed for EMIM + to enter the micropore at low surface charge. These effects are not dominant in the mesopores where the more rigid and effectively smaller TMA + ions are able to enter the pore and contribute to the increase of the capacitance.
Mixed ionic liquid electrolytes have not been widely studied, although it has been shown that they can increase the capacitance compared to the neat IL. The complex interactions arising between different molecules in the mixed IL electrolyte can change the charge densities inside the pores and thus affect the capacitance. Here, we have shown the microscopic origin for this change in the mixed IL through a detailed analysis of our MD simulations, which has not been previously done, to the best of the authors' knowledge. Therefore, we hope that our MD simulations can provide a valuable insight into mixed IL supercapacitors and can help to design new electrolytes and devises with improved properties.

Data Availability
All data needed to evaluate the conclusions in the paper are present in the paper and the Supplementary Materials. Additional data related to this paper may be requested from the corresponding authors.

Disclosure
Aleksandar Y. Mehandzhiyski's current address is Laboratory of Organic Electronics, ITN, Linköping University, 60174 Norrköping, Sweden. Candy Anquetil-Deck's current address is Norwegian University of Science and Technology, Department of Energy and Process Engineering, Høgskoleringen 5, N-7491 Trondheim, Norway.