Kinetics of Ion Transport in Ionic Liquids: Two Dynamical Diffusion States

Using classical molecular dynamics simulations, we investigate the mobility of ions in [Bmim][TFSI], a typical room temperature ionic liquid. Analyzing the trajectories of individual cations and anions, we estimate the time that ions spend in bound, clustered states, and when the ions move quasi-freely. Using this information, we evaluate the average portion of free ions that dominate conductivity. The amount of thus defined free ions comprises 15-25%, monotonically increasing with temperature in the range of 300-600 K, with the rest of the ions being temporarily bound, moving rather in local potentials. The conductivities as a function of temperature, calculated from electric current autocorrelation functions, reproduce reported experimental data well. Interestingly, for free ions the Nernst-Einstein relationship between the mobility and diffusion coefficient holds fairly well. In analogy with electronic semiconductors, one can speak about an ionic semiconductor model for ionic liquids with valence (or excitonic) and conduction band states for ions, separated by an energy gap. The obtained band gap for the ionic liquid is, however, very small, about 0.026 eV, allowing for easy interchanges between the two dynamic states.

The rediscovery of room temperature ionic liquids (RTILs)solvent free electrolytes not freezing at room temperaturewas a revolution in chemistry (1). RTILs have been appraised first not as electrolytes but as almost universal solvents for synthesis and catalysis (1), dissolving depending on their composition, practically any substance of interest. In the past we could 'count' solvents, but now with RTILs there are practically an unlimited number of them. Easily mixable with each other, 'cocktails' of RTILs offered additional variety of 'designer solvents' (2,3). However, RTILs are of great interest not only to synthetic chemists, but also to electrochemists and engineers (4)(5)(6). These electrolytes sustain higher applied voltages than aqueous electrolytes, with their ions not taking part in electrochemical reactions at electrodes when the solutes dissolved in them already do, which can speed up electrocatalytic reactions of the solutes. Together with nonvolatility and their sustainability to higher temperatures, RTILs were particularly interesting as electrolytes for electroctatalysis as well as in various energy-related devices, from supercapacitors to fuel cells and even batteries (4)(5)(6). For a physicist, RTILs comprise a unique kind of matter: room temperature dense ionic plasmas with strong Coulomb correlations.
Understanding RTILs at electrified interfaces, necessary for any electrochemical applications, required reconsideration of the theory of the electrical double layer (7). The first major attempt to do this, with discussions of avenues of further development, was presented in Ref. (8). Based on the simplest possible mean-field theory, it was demonstrated that RTILs could have dramatically different voltage dependence of the double layer capacitance from that prescribed by Gouy-Chapman theory of dilute electrolytes (9,10). 'Bell' and 'camel' shaped capacitance-potential curves, instead of Gouy-Chapman U-shape, were obtained because of a restriction on the local increase of ionic concentration in response to electrode charging. The key parameter there (and practically the only one) was the ratio of average concentration of ions in the bulk to the maximum possible local concentration of ions, γ , referred to as compacity (7). The smaller γ, the stronger the increase of the concentration of counterions in the double layer could be.
At γ =1/3 the theory predicted a crossover from the bell to a camel shape, in which the capacitance first grows and then decreases as the concentration of ions saturates at the level of maximal possible concentration.
Ref. (8) did not dwell much on the specific value of γ, but rather considered it to be a parameter of the liquid, which could be related with the free volume due to voids between ions in such complex fluids. Skipping the history of the following development of the theory of the double layer in RTILs (7), we consider the latest views on the meaning of γ.
These were influenced by two factors. Many estimates have shown that the volume portion of natural voids in RTILs is pretty small, just a few percentage (11,12). Had the value of compacity been due to the physical voids between ions, γ would not be less than 0.95. The next and even stronger impetus for reconsidering the meaning of γ came from a 'disruptive' publication (13). In that work, using surface force apparatus (SFA) the Israelachvili group measured forces between a spontaneously charged mica surface and gold with an applied voltage separated by RTILs. A striking result was reported: at large distances between the surfaces the force seemed to obey the law of an exponentially decaying repulsion of the two overlapping electrical double layers, but with an extraordinary long decay length earlier observed only for extremely diluted electrolytes. A conjecture was raised to rationalize the data (14), that pure RTILs are effectively dilute electrolytes, in which most of the ions are bound into neutral clusters (ion pairs or larger aggregates), with only a minute amount of ions participating in screening. Once adopting this picture, the rest is straightforward: (i) Large clusters of ions are majorly neutral, and have internal polarizability that provides dielectric screening to the 'free' ions; (ii) As, hypothetically, only free ions can participate in screening and contribute to the formation of electrical double layers, a very small number of such ions would give rise to extra-long screening length.
Initially this conjecture met resistance. Later Perkin et al. obtained similar results measuring structural forces between mica surfaces in contact with simple electrolytes and ionic liquid solutions of different concentrations (15). Israelachvili, Perkin, Atkin, Rutland, and co-workers have then published a mini-review about this phenomenon (16). Perkin's group has put these findings in a broader picture of nonmonotonic dependence of the screening length on the electrolyte concentration in concentrated electrolytes (17). Indeed, in dilute ionic solutions most ions are free and contribute to screening, thus the increase of their concentration makes the Debye length shorter, but with further increase of concentration more ions get bound into pairs and larger clusters and thus stop contributing to Debye screening, the latter increasing again. Such ion pairing/clustering concepts in electrolytic solutions were developed from old to modern times (18)(19)(20).
But is it true that most of the ions are indeed clustered in pure RTILs? Various estimates (19,21,22) suggested that the expected degree of clustering is not enough to explain the observations of Gebbie et al. in such a simple way. Thus, the fact that both Perkin's and Atkin & Rutland's groups have approved the ultra-long-range forces between charged surfaces, still may not exclude that the explanation of these observations lies somewhere else. A number of groups worldwide are involved in a contest to find such explanation. But before looking for alternative explanations, it would be helpful to identify whether we can indeed speak about the two states of ions, what is the exchange/balance between them, and how they contribute to conductance of RTILs, as well as to the formation of the double layer.
The cluster concept in RTILs is actually not new. In the paper by Hu and Margulis (23) the idea of long living clusters of ions has been exploited for the interpretation of the so called "red-edge effect" (REE) observed in the study of fluorescence of the organic probe 2-amino-7-nitrofluorene dissolved in 1-butyl-3-methylimidazolium hexafluorophosphate ([Bmim][PF 6 ]). The focus of that work was on establishing the relation between REE and dynamic heterogeneity in RTILs, the latter shown to be crucial for the interpretation of the data. But in the course of that study, interesting features were revealed concerning the modes of ion transport. Analyzing the van Hove self-diffusion correlation function and comparing them with those predicted by classical Fickian (Gaussian) diffusion (24), they found that most ions diffuse much slower than expected from the Gaussian diffusion, but there are a group of ions that diffuse much faster. It was also concluded that there is a poor correlation between motions of the two kinds of ions, whereas in each group the correlations are strong.
Next, in some protic RTILs, using electrospray ionization mass spectrometry (ESI-MS), Kennedy and Drummond (25) observed the formation of aggregates of ions, the size of which depends on the nature of the cation and anion. They concluded that RTILs with strong tendency for ion clustering can be classified as "poor ionic liquids", in analogy with weak, poorly dissociating electrolytes. Generally, ESI-MS has been the principal technique employed to corroborate the ion cluster model, in which the bulk structure is depicted as a sea of polydisperse aggregates (26,27). Thus, there is sufficient experimental evidence in favor of existence of clustering in RTILs but there is no unified view on the scale of this effect, which may vary from liquid to liquid.
Independent of data interpretation, the papers of Gebbie et al. had an immediate effect.
The interpretation of compacity was reconsidered and the modified mean-field theory of electrical double layer developed (22,28). Indeed, for given environmental condition (temperature, pressure, air humidity), γ can be considered as a set parameter for each RTIL in the bulk. However, near the interface, across the electrical double layer, concentration of charge carriers may vary. Indeed, in the inhomogeneous electric field the balance from the clustered state to the free state of ions shifts in favor of the free state. Thus the distribution of potential and ions across the double layer must be considered self-consistently (22). Generally, the crossover between decaying oscillatory profiles and monotonic decay of charge density correlation functions and screened potentials with increase of concentration of ionsknown as Kirkwood crossoverhas a long history. Kirkwood first predicted it, based on his analysis of the potential of mean-force in electrolytes (31). For the primitive model of electrolyte, using Ornstein-Zernike equation, Attard has studied the analytical laws for asymptotic behavior of the direct correlation function, as well as numerically its full behavior, using hypernetted chain approximation (32). He found, within the validity of approximations used, that (i) monotonically decaying screening length is slightly increased with respect to Debye length, and (ii) at high concentrations the oscillations with the period shorter than the monotonic decay length emerge. Earlier, Kjellander and Mitchell came to similar conclusions based on their charge renormalization analysis (33). Somehow further studies should establish a detailed picture of ionic states in the bulk.
The concept of ion paring and clustering in conductivity measurements and simulations of RTILs has its own history. The conductivity roughly obeys the Nernst-Einstein relationship with the measured diffusion coefficient, but with some deviations that lead to smaller conductivities (37). The concept of ionicity suggested that this deviation is due to ion correlations and ion pairing. This interpretation has been subject to harsh scrutiny, with some authors suggesting that the concept of ion pairing and clustering in RTILs is not required to explain the deviation from the Nernst-Einstein (38,39). On the other hand, some studies were in support of the concept of ion pairing and clustering in RTILs (40).
Although there have been investigations into the properties of ion pairs in RTILs (41), there have been no reports, at least as far as the authors are aware, that estimate the fraction of ions in a state that contribute to conductivity. This is the gap we aim to bridge in the present study.
The main outstanding questions here are: (i) What is the mechanism of ionic diffusion and conduction in RTILs?
(ii) What is the portion of ions that can 'freely' move and dominate the conductance and ionic screening?
(iii) Is there contribution to conductance of slowly diffusing ions corresponding to a 'bound' state?
(iv) Is the idealized picture of two states for ions physically justified for interpretation of experimental data?
(v) How easily does the exchange between the two states occur and what is the 'band gap' between them?
The present paper is focused on answering these questions. Our method of choice is molecular dynamics (MD) simulation. It allows to trace the trajectory of each individual ion, as well as to obtain the statistical characteristics of this motion such as self-diffusion coefficients, conductivity and alike, and to investigate their temperature dependence. While in the main text, we present the results only for one example RTIL,

1-n-butyl-3-methylimidazolium bis(trifluoromethanesulfonyl)imide ([Bmim][TFSI]), in the
Supporting Information (SI) we show very similar results for different RTIL, as well as test the sensitivity of results to variation of the force-fields. The conclusions remain unchanged.
Still, the answers that we will give cannot be considered final, but they set a framework for future experiments necessary for obtaining the ultimate answers.

Results and Analysis
Trajectory density approach. MD simulations of RTIL [Bmim] [TFSI] were performed at different temperatures in bulk system (Fig. 1, simulation details could be found in Methods). Respectively in Fig. 2A and 2B the trajectory of the centre of one cation and anion, which was arbitrarily chosen, is displayed. Considering diffusion of individual ions in a crowded environment can provide decisive information on the structure and dynamics of the entire system (42). It can be seen that an ion moves for a while re-entrantly (in anticorrelated fashion) in a confined volume and then 'speeds' upundergoes correlated, persistent motion outside that volume to later get trapped in a new domain. Such picture can be viewed by denser trajectory clouds connected by sparse ones. The box counting method (43) was adopted to obtain the local trajectory density (i.e., the number of trajectory points in the cubicle divided by its volume, Fig. 2C-D), by dividing the simulation cell into elementary cubicles with a size of 0.3 nm (cf. Fig. S1). To differentiate the ion states, we conventionally assume the ion to be in a 'free state' when the local trajectory-density is smaller than the average density of all trajectories; otherwise the ion is considered to be in a 'bound state'. Therefore, the percentage of thus defined free ions, which if we neglect free voids it could equalize to γ, was computed based on the times ions spend in the free and in the bound state (having averaged the trajectories separately over all cations and all anions in the simulation). The percentage of free ions is shown to increase with temperature (Fig. 2E).
That temperature dependence is treated in terms of a Fermi-like distribution, an approximate expression for the fraction of free ions was derived in Ref. (28): where E g is the energy gap between the free and bound states, and ∆S is the entropy change when an ion transfers from bound to free state. Rewriting this equation as  kJ/(mol •K) for both cation and anion. Interestingly, it is negative, indicating that the ions have more degrees of freedom in the bound state than in the free state. That entropic change comprises relatively large contribution to the free energy of transfer: at 300 K it is 3.6 kJ/mol, and at 600 K is 7.2 kJ/mol for a cation. Correspondingly, this results in a free energy of transfer 6.1 kJ/mol (0.063 eV) and 9.7 kJ/mol (0.1 eV) at 300 K and 600 K, respectively. Within the accuracy of error bars, the intercept can vary within some 10%, which will affect these numbers rather than the trends. Similar phenomenon occurs in another RTIL [Bmim][PF 6 ] (Fig. S2).
To characterize the ion dynamics in free and bound states, the velocity autocorrelation functions (VACFs) were computed for free and bound ions in

Conductivity of RTILs. The electric current autocorrelation function (ECACF) was
calculated to obtain the RTILs conductivity, since ECACF takes into account not only the ionic motions but also the contributions of cross-correlations among ions (45). It can be found that ECACFs decay very quickly to zero within 1 ps (Fig. 4A) and the conductivity spectrum shows a higher peak shifting toward low frequency as the temperature increases ( Fig. 4B). Data for all studied temperatures can be seen in Fig. S4.
In view of the Green-Kubo integral of ECACF (44), including the contributions of ionic correlations, the electrical conductivity, Λ, is computed for the RTIL [Bmim] [TFSI] and exhibits good agreement with experimental measurements (49), as shown in Fig. 5 (the red squares vs the blue dots). Several research groups reported that the electrical conductivity would be overestimated by Nernst-Einstein equation, because it neglects ionic association and assumes that each ion contributes one elementary charge (37,44,45,50).
Such overestimation would be clearly seen in Fig. 5, comparing the conductivity obtained by full ECACFs with that calculated from Nernst-Einstein equation (the red squares vs the black circles).

Fig. 5. The electrical conductivity, Λ, of [Bmim][TFSI] as a function of temperature.
The red square shows the values calculated from the MD-simulated electric current autocorrelation function (Fig. 4A). The blue dots represent experimental data from Ref. (49). The black circles display the values computed via original Nernst-Einstein equation with diffusion coefficients shown in Figure 3B. The magenta stars are based on Eq. (2), with free ion percentage evaluated by the trajectory density method. The green diamonds are based on Eq. (2), with free ion percentage evaluated by ion pairing method.
With the above analysis of free and bound ions, let us assume that the free ions dominantly contribute electrical conductivity, and the Nernst-Einstein equation could be

Discussion and Conclusions
Mechanisms of Ion Transport in RTILs. There were different ideas about the mechanisms of ion transport in RTILs. One of the scenarios was proposed by Abbot (52,53) who assumed that in such a concentrated ionic systems, only a hole mechanism of transport would be feasible. The findings presented above cannot directly prove or disprove that hypothesis, as looking at any individual trajectory, we cannot tell how each elementary step of ion transfer proceeds. Another, often discussed, such scenario is a quasi-Grotthuss (relay) mechanism, in which the ions in each elementary act do not move far, but each ion just slightly shifting, kicking another ion to continue. Whereas this is not excluded for ions inside the clusters, our free ions can seemingly individually move considerable distances, as we can see it from their individual trajectories.
Ionic Semiconductor Analogy. As typical in condensed matter physics, in complex cases people indulge into the language of 'collective' phenomena. It is, thus, tempting to treat RTILs as an 'ionic semiconductor', in which in the ground state ions populate the ionic 'valence' (or excitonic') band, but they can be excited into a 'conductance' band (54).
However, the band gap energy that we obtained in the fits of More experiments are needed to verify this concept. Measurements of the Hall effect come to mind first. However, the latter is difficult from experimental point of view (58,59), as well as is not straightforward in interpretation (60,61).
Our last note is on possible implications for the double layer theory. If the exchange between the conduction and valence bands is much faster than the RC-time of charging the double layers, which is practically always the case, it would be legitimate to speak about the average concentration of mobile charge carriers that contribute to the formation of the electrical double layer. But if even only 20% of ions were on average free and giving the major contribution to the double layer, it is already a very concentrated ionic system, and the spatial structure of the double layer, if described in all its complexity, will presumably contain oscillating ('overscreening') part and the exponential tail ('underscreening' part).
In fact, a recent study of the temperature dependence of the differential capacitance found the shape of the capacitance curve to change as if ion pairs and clusters were breaking up with increasing temperature (28), which is in accordance with the qualitative results here (although the reported values of  therein were slightly different). Bound ions will contribute not only to effective dielectric constant of the system, but also to the spatial potential and charge distributions.

Methods
As shown in Fig. 1 Supplementary Information Fig. S1 shows the size distribution of different clusters of ion trajectory cloud ( Fig. 2A-B) in bound state. One can see that the most trajectory cloud size is around 0.3 nm. Therefore, we adopted 0.3 nm as the optimal cubicle size for the box counting method.
trajectory cloud size (nm)    As Fig. S5 shows, modified Nernst-Einstein equation, taking into account the introduced parameter of free ion percentage (see Fig. S2A), gives the electrical conductivity, in line with the approach of the Green-Kubo integral based on full ECACFs as well as the experimental data (1).