Superionic Liquids in Conducting Nanoslits: Insights from Theory and Simulations.

Mapping the theory of charging supercapacitors with nanostructured electrodes on known lattice models of statistical physics is an interesting task, aimed at revealing generic features of capacitive energy storage in such systems. The main advantage of this approach is the possibility to obtain analytical solutions that allow new physical insights to be more easily developed. But how general the predictions of such theories could be? How sensitive are they to the choice of the lattice? Herein, we address these questions in relation to our previous description of such systems using the Bethe-lattice approach and Monte Carlo simulations. Remarkably, we find a surprisingly good agreement between the analytical theory and simulations. In addition, we reveal a striking correlation between the ability to store energy and ion ordering inside a pore, suggesting that such ordering can be beneficial for energy storage.

Introduction

Confined ionic liquids show an exciting physics and play an increasingly important role in technology, finding their applications in various electrochemical devices, such as capacitive energy storage[1−4] and deionization systems,[5−7] heat-to-energy converters,[8−11] etc. Electrochemical capacitors, or supercapacitors, for instance, store energy through charge separation between their cathode and anode.[12] Since the amount of the stored charge scales with the contact area between the electrode and an electrolyte, a method to boost the energy storage is to use electrodes with large volume-filling surface area, which are provided by highly porous electrodes, particularly those containing subnanometer pores. Such pores, which can admit just a single layer or row of ions, deliver the highest achievable capacitance,[13−15] allegedly due to a superionic state emerging in a narrow conducting confinement.[16] In the superionic state, the electrostatic interactions between the ions are strongly screened,[16−19] allowing a tighter packing of counterions and easier unbinding of “ion pairs”, thus enhancing the capacitance.[16,20]

The simplest model for charge storage in subnanometer pores is a classical antiferromagnetic spin-1 Ising model with nearest-neighbor interactions in external field.[21] In this model, cations and anions correspond to ±1 spins and the external magnetic field corresponds to the potential drop between the electrode and bulk electrolyte. In one dimension (1D), the Ising model mimics a single-file pore and obeys the well-known exact solution,[22] allowing physically appealing analyses of the charging behavior.[21] In deviance of its simplicity, this model turns out to capture the essential physics of charging nanoscale pores quite well.[23,24] More complex 1D models have also been used in various contexts.[25−31]

Mapping the confined ions onto 1D lattice spin models has an advantage that such models are often exactly solvable in external field, which corresponds to the cell voltage in electrochemical devices. Such quasi-1D pores are typical for electrodes based on open carbon nanotubes (CNTs) or on a forest of closed CNTs. For the state-of-the-art graphene-based[32−37] and MXene[38−40] electrodes, slit pores are a better approximation to the reality. New physics emerges in such pores and particularly exciting is the possibility to observe phase transitions, which cannot occur in 1D systems.[41,42]

Two-dimensional lattice models are widely used in the literature to study confined or adsorbed ions, but they are mainly handled by simulations.[43−46] Recently, a three-component model has been proposed[47,48] for ionic liquids in narrow conducting slits, which can be mapped onto the Blume–Capel (BC) model, well known in the theory of magnetism.[49−53] Similarly to the Ising model, the BC model has been solved exactly only in 1D. In refs (47, 48), this model was treated with a Bethe-lattice approximation, which allowed analytical insights into the dependence of charging on applied voltage and other parameters. A natural question is, however, how accurate and reliable these analytical results are. In a broader context, the Bethe-lattice approximation is a convenient analytical tool that is frequently used to tackle various models of statistical physics, and hence the above question is of a general relevance.

Herein, we assess the accuracy of the Bethe-lattice approximation with Monte Carlo (MC) simulations of the model of ref (47) on the square, honeycomb and, in a few cases, triangular lattices. In addition, with the Bethe-lattice calculations, we discuss which of these three lattices is more likely to be realized in off-lattice simulations and experiments. With MC simulations, we gain additional insights into the structural transformations of ionic liquids in narrow slit confinements and investigate how ordering of ions is related to energy storage.

Of course, our model does not account for many properties of real electrodes and electrolytes, such as the carbonic nature of pore walls, dispersion interactions, ion and solvent polarizability, etc. (see Section II.C in ref (48) for a detailed discussion of model limitations). In this sense, the model is probably too idealized to directly compare with experiments using the currently available microporous carbons, although it might be possible to do so with future nanostructured electrodes. Nevertheless, we expect some qualitative features of the charging behavior and stored energy–ion ordering relations, revealed in this article, to be generic and model-independent.

Model and Methods

Model

To study charge storage in ultranarrow slit nanopores that are so narrow that they admit only one layer of ions, we consider a model Hamiltonian defined on a two-dimensional lattice[47,48] (Figure )Here, n+ = 1 (n– = 1) means that site i is occupied by a cation (an anion) and n+ = n– = 0 means that site i is vacant or occupied by solvent (configuration n+ = n– = 1 is prohibited due to hard-core exclusion); ⟨ij⟩ denotes nearest-neighbor sites and I > 0 is the interaction energy between two neighboring ions. For ions confined between two metallic plates I ≈ ϕ (a), where a is the lattice constant (≈the ion diameter) and[16]with K0(x) being the modified Bessel function of the second kind of order zero, L the slit width, e the elementary charge, ε a dielectric constant inside a pore, which comes from electronic and rotational degrees of freedom of ions and due to a solvent, if present. The value of ε is not precisely known, but we expect it to vary from 2 for simple ions to 5 or more for more complex, bulky ions or in the presence of a solvent.

Figure 1

Model for ions in a narrow slit nanopore. (a) Ions are confined into a narrow metallic slit. Ion diameter a = a± is the same for cations and anions. Potential difference u is applied at the pore walls with respect to the bulk electrolyte (not shown). The energy of transfer of ions from the pore into the bulk w = w± is the same for both ions (see text). (b) Arrangement of ions on square, honeycomb, and triangular lattices, used in Monte Carlo (MC) simulations of model (1). Blue and orange spheres represent ions, and gray spheres represent void/solvent. Interaction energy between the same type of ions is I, and that between oppositely charged ions is −I (see eq ). (c) Fragment of the Cayley tree with coordination number q = 4, corresponding to the square lattice. (d) Ordered phases on the square lattice. Ordered phase I consists of an equal amount of cations and anions occupying two sublattices. In ordered phase II, the counterions occupy one sublattice while the other sublattice is occupied by solvent/void.

Equation gives coupling constants I that vary from just a kBT for large ions and narrow slits (L/a ≈ 1) to a few kBT for small ions or wide slits. I increases for increasing L at constant a and decreases for increasing a at constant L (Fig. S1). Note that we consider 1 < L/a < 2, that is, the confined ions form a single layer. For L/a values close to but smaller than 2, the next-nearest interactions may become important.[48] It will be interesting to study such effects in future work.

In eq , “external” fields h± arewhere u is the applied potential (with respect to the bulk electrolyte outside of the pore) and w± is the energy of transfer of a ± ion from an empty into the bulk, which includes the image forces and other interactions of the ions with the pore walls.[16] We note that the definition used here differs by sign from the re-solvation energy used in other works.[16,54,55] We assume w+ = w– = w and note that the results for an asymmetric ionic liquid (w+ ≠ w–) can be obtained by shifting the applied voltage by −(w+ – w–)/2e and taking the transfer energy equal to (w+ + w–)/2.

Based on quantum-mechanical calculations and molecular dynamics simulations,[56] Lee et al.[57] have estimated that for slit pores, the nonelectrostatic contribution wnon to w ranges from −5 to 45 kBT. The ion’s self-energy wself (due to image charges[16]) is negative and varies from a few negative kBT’s for wide slits to about −50 kBT for narrow slits (Figure S2). Thus, the total ion’s transfer energy, w = wself + wnon, can span a wide range of values from negative to positive around zero. We note that changing the slit width or temperature affects both βI and βw.

Monte Carlo Simulations

We performed two-dimensional (2D) Monte Carlo (MC) simulations of model (1) on square, honeycomb, and triangular lattices (Figure ). Periodic boundary conditions were applied in both directions. The lattice size was 32 × 32 sites in all simulations, but in selected cases, we carried out simulations for smaller and larger lattices to study the finite-size effects. We observed significant deviations only in the vicinity of second-order transitions, as one might expect (Figures S4, S6, and S8). The system was equilibrated with 2 × 105 MC steps, each consisting of M = 322 single-spin updates, and the averaging was performed over 8 × 105 MC steps.

Bethe-Lattice Solution

The free energy of the system iswhere the sum runs over all possible configurations of the occupation numbers n+, n– on a given lattice. Within the Bethe-lattice approach, eq is evaluated on a Bethe lattice with the same coordination number q as of the original lattice (Figure c). The partition function on the Bethe lattice was computed exactly. For details of the calculations, see ref (47) (nonpolarized electrodes) and ref (48) (voltage dependence, q = 3). In this work, we took q = 4, 3, and 6, corresponding to square, honeycomb, and triangular lattices used in the MC simulations.

Ordered and Disordered Phases

Recently, ref (48) has demonstrated a rich phase behavior of model (1) (with q = 3 neighbors), involving direct and re-entrant symmetry-breaking phase transitions between “ordered” and “disordered” phases. The disordered phase is a homogeneous mixture of ions and voids, while ordered phase means that the ions of one type predominantly occupy one of the two sublattices. We will see that at low applied voltages, the ordered phase consists of an equal amount of cations and anions (ordered phase I) so that the pore is neutral or only weakly charged. In the ordered phase at high potentials, one sublattice is occupied by counterions and the other sublattice is mostly free of ions (ordered phase II). Typical ion configurations in these phases are schematically shown in Figure d.

Thermodynamic Characteristics and Charging

Quantities directly accessible from the Bethe-lattice calculations and MC simulations are the average number of ions ⟨n±⟩ in the pore. It is convenient to discuss charging in terms of the total ion density and the accumulated chargeNote that by definition, these quantities are measured per lattice site and to transform them to the corresponding quantities per surface area, one has to divide them by the area per site, which is a2, 3a2√3/4, and a2√3/2 for the square, honeycomb, and triangular lattices, respectively.

Having calculated Q(u), one can compute the differential capacitanceIn MC simulations, however, it is more convenient to evaluate the capacitance from charge fluctuations[58]Having the capacitance, one can calculate the energy stored in a nanopore at potential difference uSimilar to charge and density, the capacitance and stored energy are measured per lattice site.

We shall also calculate the charging parameter[59,60]which describes charging mechanisms taking place at applied voltage u. It is not difficult to see that XD can be expressed in terms of charge-density fluctuationsWe have used eq to compute XD in MC simulations.

Finally, we shall study the behavior of density fluctuations, related to the isothermal compressibilityThis quantity describes how liable the system is to compression during charging, which may be relevant to electroactuators[61] and studies of electrode deformations.[62]

Results and Discussion

Phase Diagram for Nonpolarized Electrodes

Figure shows a phase diagram spanned in the plane of ion–ion interaction energy I and ion transfer energy w for a nonpolarized electrode (i.e., potential difference u = 0 with respect to the bulk electrolyte). The phase diagram consists of disordered and ordered phases separated by the lines of first-order and second-order phase transitions, which meet at a tricritical point (filled circle in Figure ). As discussed, the disordered phase is a homogeneous mixture of ions and solvent/voids, whereas the ordered phase (type I) consists of an equal amount of cations and anions occupying two sublattices of the square lattice (Figures d and S3). For the phase diagram on the honeycomb lattice, see Figure S5.

Figure 2

Phase diagram of a superionic liquid in a nonpolarized ultranarrow slit. The results of the Bethe-lattice calculations with coordination number q = 4 are shown by lines, and the square symbols denote the results of Monte Carlo simulations on a square lattice. The filled circle shows the tricritical point obtained within the Bethe-lattice approach. The black solid line and filled squares denote second-order continuous transitions, and the black dashed line and open squares show first-order discontinuous transitions between a disordered phase and the ordered phase I. For schematics of ordered and disordered phases, see Figure . The thin vertical lines show the values of the transfer energy chosen to study the charging behavior (Figures and 4). For the phase diagram on the honeycomb lattice, see Figure S5.

Figure 3

Phase behavior and charging of pores with transfer energy w = 0. (a) Phase diagram in the plane of applied voltage u and ion–ion interaction energy I. The solid line denotes a line of second-order phase transitions separating the disordered phase and the ordered phase (type I). The thin vertical lines show the values of βI used in the remaining panels. The diagram has been determined using the Bethe-lattice approach. (b) Charge on sublattices A and B as a function of voltage. The charge on both sublattices is the same for βI = 0.4. For βI = 2, there is an ordered phase for βue ≲ 8.5, in which the charges on two sublattices have opposite sign but the same magnitude. (c) Total ion density, (d) capacitance, (e) stored energy, and (f) compressibility as functions of applied voltage. The inset in (f) shows the charging parameter XD (Equation and 8b). The lines are the Bethe-lattice results, and the symbols denote the results of MC simulations. (g) Snapshots from Monte Carlo simulations with βI = 2. For finite-size effects, see Figure S6, and for the results on the honeycomb lattice (q = 3), see Figure S7.

Figure 4

Phase behavior and charging of pores with transfer energy βw = −4. (a) Phase diagram in the plane of applied voltage u and ion–ion interaction energy I. The solid and dash lines denote the lines of second-order and first-order phase transitions, respectively, separating the ordered and disordered phases. The thin vertical lines show the values of βI used in the remaining panels. The diagram has been obtained using the Bethe-lattice calculations. (b) Charge on sublattices A and B, (c) total in-pore ion density, (d) capacitance, (e) stored energy, and (f) compressibility as functions of voltage. The inset in (f) shows the charging parameter XD (Equation and 8b). The lines are the Bethe-lattice results, and the symbols denote the results of MC simulations. (g, h) Snapshot from Monte Carlo simulations for βI = 2 and 3. For finite-size effects, see Figure S8, and for the results on the honeycomb lattice (q = 3), see Figure S9.

From MC simulations, the locations of the transitions were determined by analyzing the density profiles on two sublattices (Figure S3). While for the discontinuous transitions far from the tricritical point, this procedure leads to an excellent agreement with the Bethe-lattice calculations, for the second-order transitions, there are some deviations. In particular, in the MC simulations, the transitions are shifted toward higher values of βI (Figures , S3, and S4a). This is probably because the Bethe-lattice approach underestimates the effect of fluctuations in the ordered phase. Nevertheless, the agreement between the MC simulations and analytical Bethe-lattice calculations is remarkable.

The diagram in Figure shows the states of systems characterized by different values of coupling constant I and transfer energy w. We now select a few typical systems from this diagram and investigate how they respond to a voltage applied to the pore with respect to the bulk electrolyte.

Charging

Using the Bethe-lattice approach, ref (48) has demonstrated a complex charging behavior of model (1) for coordination number q = 3, corresponding to the honeycomb lattice. This behavior is essentially determined by the ion transfer energy w and can be divided into four zones related to the number of phase transitions that the system undergoes during charging. Here, with the Bethe-lattice approach, we find a qualitatively similar behavior for q = 4. We assess these Bethe-lattice results with MC simulations on the square lattice, which allow us to gain additional physical insights into the transformation of ionic structure under an applied voltage. These results are discussed below. The corresponding comparison for the honeycomb lattice is presented in the Supporting Information (Figures S7 and S9).

Positive Transfer Energy of Ions

Positive transfer energies correspond to ionophilic pores that are occupied by ions at zero potential difference u. The charging behavior depends on whether the ions are in the disordered or ordered phase at u = 0. In the case of the disordered phase (βI = 0.4 in Figure ), the charging proceeds continuously and the ions remain in the disordered phase at any applied voltage (vertical black dash line in Figure a). The ion density increases with u due to counterion adsorption, and the capacitance exhibits a small peak close to u = 0 before decaying to zero (black lines in Figure c,d). Consequently, the stored energy increases and saturates at high voltages (Figure e). The charging mechanism consists of a combination of swapping and electrosorption (XD > 0, the black line in the inset of Figure f). Interestingly, the compressibility has a peak at u = 0 but decreases quickly to zero, similarly to the capacitance (Figure f). At all potential differences, we found an excellent agreement between the Bethe-lattice approach and the MC simulations (black lines and squares in Figure ).

If the system at u = 0 is in the ordered phase, it undergoes a second-order phase transition to the disordered phase upon increasing the applied voltage (βI = 2 in Figure ). In the ordered phase type I, two sublattices (A and B) are occupied by an equal amount of cations and anions (Figure b,g). The ions remain in the ordered state at low voltages so that the accumulated charge and hence the capacitance and stored energy practically vanish in this regime (Figure d,e). Upon further increase of voltage, the charging commences by expelling the co-ions from the pore, which implies that the total ion density decreases and the charging parameter XD < 0 (Figure c and the inset in Figure f). Both the capacitance and the compressibility show a peak at the transition (Figure d,f). In the MC simulations, the locations of the transitions are slightly shifted toward the ordered phase, compared to the Bethe-lattice predictions. This is likely because the Bethe-lattice approach underestimates the effect of fluctuations and hence overestimates the stability region of the ordered phase.

Interestingly, the stored energy at high applied potentials (viz., in saturation) is about 4 times higher for βI = 2 than for βI = 0.4 (Figure e). This is because for βI = 2, the system is in the ordered state at u = 0, which persists until a sufficiently high voltage is applied. Thus, the charging is effectively shifted to higher voltages, which leads to higher stored energies. This is similar to charging ionophobic pores, which can provide much higher stored energies compared with ionophilic pores as long as the pore ionophobicity shifts the charging process to higher voltages.[63]

Negative Transfer Energy of Ions

A negative transfer energy means that an ion prefers to stay outside of a pore. However, this does not necessarily imply that the pore is empty at zero voltage, as the pore occupancy depends also on the in-pore ion–ion interaction energy βI. A high βI can compensate the unfavorable transfer energy, thus promoting the ion pairs to go inside the pore. Increasing βI therefore increases the pore occupancy.

The case of low βI has been studied in ref (48) with the Bethe-lattice approach (for coordination number q = 3, see Figure S9). Here, we found a qualitatively similar behavior in the case of a square lattice (q = 4; see βI = 1, 2 in Figure ). For βI = 1, the charging proceeds continuously with u. For βI = 2, there are two phase transitions that the system undergoes upon increasing voltage: first from the disordered phase to the ordered one, and then back to the disordered phase. The ordered phase is type II and consists of counterions occupying one sublattice, whereas the other sublattice is free of ions (βI = 2 in Figure b,g). Interestingly, a similar sequence of the voltage-induced transitions has been reported for an adsorbed layer of butylmethylimidazolium-hexafluorophosphate (BMIM-PF6) ionic liquid, forming fluid-like (disordered) and lattice-like (ordered) phases at a graphite electrode.[64]

We found that the MC and Bethe-lattice approaches are again in good agreement, except for the close vicinity of the transitions. For βI = 2, the simulations show that the stability of the ordered phase is reduced, compared with the Bethe-lattice results, which is likely due to fluctuations that are not fully accounted for in the Bethe-lattice approach.

The case of high βI’s has been omitted in ref (48), but it is interesting in that the type of ordering changes in the course of charging. At zero voltage, the pore is filled with ions, which form type I ordered phase (βI = 3 in Figure ). With increasing the voltage, the counterions remain on their sublattice, while the co-ions gradually leave the pore. Perhaps surprisingly, this process proceeds without a transition. Remarkably, a similar behavior has been observed in molecular dynamics simulations of charging ultranarrow slit nanopores, where the charging proceeded via “melting” of an interface between the two types of “quasi-ordered” ionic liquid states.[65]

Similarly to the case βw = 0 (Figure ), the stored energy in saturation is much higher if the system is in the ordered state at zero voltage (βI = 3 in Figure e).

Square, Honeycomb, or Triangular?

We have discussed the charging behavior of superionic liquids using the lattice model, eq , with the ions arranged on the honeycomb and square lattices. Although both lattices give rise to a virtually identical qualitative behavior, there are quantitative differences (compare, e.g., Figures and S7). A natural question thus arises as to which of these two types of positional ordering (number of neighbors) could be realized in experimental systems or off-lattice simulations. To address this question, we used the Bethe-lattice approach to compare the free energies (eq ) of superionic liquids during charging on these two lattices. In addition, we studied charging on a triangular lattice (the number of nearest neighbors q = 6), which provides the densest packing of ions. However, the triangular lattice is tripartite, that is, it contains three sublattices and does not allow ordering of ions on two sublattices. Since ordering on the triangular lattice can only be realized at specific concentrations of ions and solvent, here, to avoid these complications, we restrict our considerations to the disordered phase where the sublattice division plays no role. We note that also for the triangular lattice we found a good agreement between the Bethe-lattice calculations and MC simulations (Figures S10 and S11).

For an ionophilic pore, the system on the triangular lattice has the lowest free energy at all voltages (Figure a). Qualitatively, the ion density and the capacitance behave similarly on all lattices (Figure b,c), but, as mentioned, the triangular lattice provides the highest density.

Figure 5

Square, honeycomb, or triangular? (a) Free energy per surface area, F; (b) 2D ion density, ρ; and (c) capacitance per surface area, C, for honeycomb, square, and triangular lattices obtained by the Bethe-lattice calculations with the coordination numbers q = 3, 4, and 6, respectively. The transfer energy of ions w = 0 and the ion–ion interaction constant βI = 0.2 correspond to an ionophilic pore, occupied by ions at zero voltage. (d–f) Same as (a–c), but for the parameters βw = −4 and βI = 0.8 corresponding to an ionophobic pore that is nearly free of ions at no applied voltage. The dashed lines denote “metastable” branches having higher free-energy densities. For clarity, only the branches with the lowest free energy are shown in (e, f). The thin vertical lines denote the transformations between the square and triangular ordering.

For an ionophobic pore, there is a re-entrant “transition” between the triangular and square lattices (Figure d). The square lattice has the lowest free energy for intermediate voltages because, for an increasing ion density, it allows the counterions to avoid the unfavorable direct contact with each other by having fewer nearest neighbors. At sufficiently high voltages, the triangular lattice provides the densest packing and hence the lowest free-energy density. It is interesting to note that hexagonal-symmetry ordering (i.e., the triangular lattice) has been reported in simulations of BMIM-PF6 at a graphite electrode.[66]

It is instructive to compare our results with recent 2D off-lattice simulations by Schmickler and Henderson, who revealed spectacular weblike patterns of an ionic liquid in the regime of weak ionophilicities.[67] These authors also reported on enhanced fluctuations, characteristic of criticality, in some range of parameters. Interestingly, they saw a tendency for ions to form a square lattice at zero voltage, while our Bethe-lattice analysis suggests that hexagonal-symmetry ordering is preferential. We note, however, that our analysis was applied in the regime of the disordered phases (i.e., the same average ion densities on both sublattices), which corresponds to very low ion–ion interaction energies (βI = 0.2). The slit width considered in ref (67) was 5a, where a is the ion diameter. For such wide slits, the coupling constant reaches the value of about 200 kBT, according to eq . Moreover, the interactions between the next-nearest and further neighbors become non-negligible (e.g., ϕ(2a) ≈ 70 kBT) and hence the results of our simple model cannot be directly compared to the case of ref (67).

Conclusions

We have studied the charging behavior of superionic liquids in a narrow slit confinement using the recently introduced three-component lattice model.[47,48] We solved this model exactly on a Bethe-lattice and carried out the corresponding Monte Carlo simulations on the square, honeycomb and triangular lattices. Surprisingly, we found a remarkably good agreement between the analytical results and MC simulations, except for the regions in a close vicinity of second-order transitions. This result is exciting and encourages the application of the analytically tractable Bethe-lattice approach to other systems of electrochemistry and statistical physics.

With MC simulations and Bethe-lattice calculations, we confirmed the rich phase behavior of superionic liquids under an applied voltage, as found in an earlier work.[48] In addition, we revealed that the energy stored in a pore is a sensitive function of the thermodynamic state. If the in-pore ions are in the ordered state at zero voltage, the system remains in this state until a sufficiently high voltage is applied to break the order, which effectively shifts the charging process to higher voltages and can lead to even a few-fold enhancement of the energy storage. This is similar to the effect of charging ionophobic pores that are free of ions in the nonpolarized state. Indeed, in both cases, the charging is shifted to higher voltages: for ionophobic pores, this is because of an energy barrier for the ions to enter the pore;[63] for ordered ionic structures inside nonpolarized pores, it is due to the difficulty to break this structure before enriching the pore interior with counterions. Creating a difficulty to do more work to store more energy is a general concept.

Thus, for charging at high voltages, the initial in-pore ordering of ionic liquids might be beneficial for energy storage. Although counterintuitive, this conclusion may appear to be general and model-independent.