Lattice model of ionic liquid confined by metal electrodes

Room Temperature Ionic Liquids (RTILs) have attracted a lot of attention due to many potential applications.¹ RTILs present high capacitance and fast charging,² substituting dielectric materials in supercapacitors.^3–6 They are stable under a large voltage window, sometimes reaching 0-6 V, which makes them ideal for energy storage. They have been investigated for possible applications in microreactors,⁷ in solar cells,^8,9 and other renewable energy devices,^10–12 making the understanding of these complex fluids of paramount importance.

Ionic liquids have inherently strong electrostatic interactions between the ions and very high packing fractions. These characteristics make them very difficult to describe using existing theoretical methods. However, these same features make these systems particularly interesting on purely fundamental grounds¹³ as they stand at the frontier of statistical mechanics¹⁴ demanding the development of more sophisticated theoretical methods.^15–21 

The most widely used theoretical approach to study ionic liquids is based on the modified Poisson-Boltzmann equation (mPB).^22–26 The mPB equation approximately accounts for the finite size of both ions and solvent and has been used to calculate differential capacitance of RTILs. However it has been shown^27–29 that the mPB equation does not accurately predict the double layer structure of electrified interfaces. Therefore more sophisticated approaches based on the Density Functional Theory (DFT)^30–34 must be used. Unfortunately, such theories are numerically very demanding and it is also very difficult to account for strong electrostatic correlations within the DFT formalism.³⁵ An alternative approach to study strongly correlated Coulomb systems is to use computer simulation. However, even in bulk, the long range nature of Coulomb interactions makes numerical simulations very demanding.^36–38 The fundamental problem is that because of the long range interaction, one cannot use standard periodic boundary conditions. To properly explore the thermodynamic limit, one needs to periodically replicate the simulation box so that the ions in the main simulation cell also interact with the ions in an infinite set of replicas. To efficiently sum over replicas, Ewald summation techniques^39–44 have been developed. Confinement of Coulomb systems presents additional difficulties,⁴⁵ requiring more specialized techniques.^46–53 The computational cost of summation over replicas increases dramatically due to appearance of special functions and slow convergence in quasi two-dimensional systems.⁵⁴ In various applications, the simulation box must be made sufficiently large to achieve a bulk-like regime in between the surfaces,²³ requiring a large number of particles, which further slows down simulations.²⁹ If the confining surfaces are polarizable, there are additional complications connected with the induced surface charge.^55–57 The usual methods for simulating metallic surfaces of electrodes are computationally very expensive, relying on a minimization procedure to calculate the induced surface charge at every simulation time step.^58–61 Alternatively, there are electrostatic layer correction methods (ELC) or methods that explicitly sum infinite series of image charges.^62–64 

Recently, we introduced a new simulation method which allows us to easily calculate sums over replicas.⁶⁵ The formalism is based on periodic Green’s functions and was applied to a lattice model of a symmetric RTIL liquid; see Fig. 1. The algorithm allows one to rapidly calculate the induced surface charge, permitting efficient simulations in a constant electrostatic potential ensemble.⁶⁶ Here we pursue further this approach to calculate the differential capacitance curves for ionic liquids both symmetric and asymmetric in ionic size and charge. The great advantage of the lattice model simulations is that ionic interactions can be precalculated at the beginning of simulations, making such simulations substantially faster.

The paper is organized as follows: First, we briefly review Green’s functions for ions confined by metal surfaces. Next, we discuss Monte Carlo (MC) simulations in a constant electrostatic potential ensemble. We then present the density profiles and the differential capacitance curves for various ionic liquids. The paper will conclude with the discussion of results and a perspective for future work.

We start by briefly reviewing the calculations of electrostatic potential produced by an ion of charge q confined by parallel infinite grounded conducting surfaces located at z = 0 and z = L; for more details, the reader is referred to Ref. 65. The electrostatic potential can be conveniently written in cylindrical coordinates,⁶⁷ 

where (0, z₀) is the coordinate of the ion, (ρ, z) is the observation point, k_n = nπ/L, and K₀ is the modified Bessel function of second type. Unfortunately, this expression has convergence difficulty when ρ → 0. In this case, we can find a different representation of the Green’s function⁶⁷ 

where J₀ is the Bessel function of order zero. This equation is well behaved when ρ → 0, as long as z ≠ z₀. The presence of an ion between the grounded conducting surfaces will result in an induced surface charge which can be calculated from the discontinuity of the normal component of electric field at the surface. Integration over the surface of each electrode allows us to obtain the induced total surface charge.⁶⁵ We find

for the left and right electrodes, respectively. If the electrodes are not grounded, but are held at a constant potential difference ψ₀, then there is an additional contribution to the electrostatic potential,

The simulations are performed in the NVT ensemble at a fixed electrostatic potential difference ψ between the electrodes.^65,66 The partition function assumes the form

where β = 1/k_BT and the surface charge on the left and right electrodes within the simulation cell is ∓Q, respectively. Since in this ensemble the surface charge on the electrodes can fluctuate, the calculation of the differential capacitance of the system is straightforward

where A = L_xL_y is the area of the electrode in the simulation cell, which has volume V = L_xL_yL. Note that in order to perform simulations in the fixed surface potential ensemble, we need to know the electrostatic energy at a fixed surface charge, E(Q). The charge distribution will not be uniform over the surface of the electrodes and will respond to the ionic motion.

The simulation cell is periodically replicated in x and y directions. The electrostatic energy⁶⁵ inside the simulation cell for a given surface charge ±Q is

with the periodic Green’s function constructed using Eq. (1)

where m’s are integers corresponding to periodic replicas of the system. U_s(r_i) is the self energy of each ion calculated using the limiting process,

which yields, with the help of Eq. (2), the following expression:⁶⁵ 

For a fixed surface charge Q, the surface potential ∓ψ₀/2 on each electrode will fluctuate as a result of ionic motion. Since the system is charge neutral, the surface potential⁶⁵ for a given surface charge and ionic distribution inside the simulation cell can be easily calculated using Eqs. (3) and (4),

Now we are in a position to perform MC simulations in the constant electrostatic potential ensemble, in which the surface charge Q fluctuates in accordance with Eq. (5).

The simulation cell has volume V = L_xL_yL, with L_x = L_y = 64 Å and L = 240 Å in the case of symmetric ionic liquids and L = 160 Å otherwise. The ionic liquid is confined in the region −L_x/2 < x < L_x/2, −L_y/2 < y < L_y/2, and 0 < z < L. The Bjerrum length, defined as λ_B = q²/ϵk_BT, assumes two values: λ_B = 7.2 Å which is appropriate for electrolytes dissolved in pure water at room temperature and λ_B = 38.4 Å, which is suitable for RTILs.^68–70 The ions are constrained to move on a lattice with lattice spacing a = 4 Å, a = 8 Å, or a = 16 Å, depending on the system studied. We start by considering symmetric ionic fluids with spherical ions of diameter equal to the lattice spacing and charge ±q, where q is the proton charge. The compacity factor γ is defined by the ratio between lattice sites occupied by the particles and the total available lattice sites in the simulation box; see Fig. 1. In order to properly sample the phase space, we use both short- and long-range movements. The differential capacitance is calculated using 4 × 10⁵ uncorrelated samples after equilibrium has been achieved. Fig. 2 shows the change in the form of the differential capacitance between strong and weak coupling regimes. Note that the differential capacitance is symmetric with respect to ψ → −ψ. In the weak coupling electrolyte regime, the differential capacitance has a characteristic camel-back shape, while in the strong coupling regime, it has a bell-shaped form. Contrary to the predictions of mPB theory,²² the transition between camel-shaped and bell-shaped regimes does not occur at the universal value γ = 1/3 but instead depends on both γ and λ_B, see Fig. 3.

For ionic liquids (strong coupling regime) with charge asymmetric 2:1 ions, Fig. 4 shows appearance of a second peak at the intermediate applied voltages, if the system is sufficiently dense. This is similar to what has been found in continuum simulations.⁶⁰ We see, however, that the height of the secondary peak does not scale with the surface area of the simulation box A = L_xL_y; therefore, there is no structural phase transition in the local ordering of ions near the electrode, contrary to the suggestion in Ref. 60. Figs. 5 and 6 show that the secondary peak is correlated with the appearance of additional structure in the ionic density profiles—the peak is absent for γ ≤ 1/10, see Fig. 5, when profiles show less layered structure, Fig. 6. In both cases, with and without second peak, γ = 4/10 and γ = 1/10, respectively, we find that the first layer overscreens the charge on the electrode—the charge on the cathode is ≈−60q, while the contact layer has charge of ≈+80q. Such strong correlational effects clearly cannot be captured by mean-field theories, such as mPB. Finally we consider the differential capacitance of a 1:1 ionic liquid with size asymmetric (two-to-one) ions, see Fig. 7. Following Ref. 71, a large ion excludes 93 vertices and a small ion 27, of the total available in the simulation box. The compacity factor is then γ = (93N₊ + 27N₋)/N_box, where N_box is the total number of lattice sites. Fig. 8 shows that size asymmetry leads to a reduction of the maximum differential capacitance. The magnitude of this reduction is similar to the one found in symmetric systems, with ions of radius R = 8 Å, Fig. 2(d). The capacitance curve of an asymmetric RTIL has also significantly more structure.

We used MC simulations in a constant electrostatic potential ensemble to calculate differential capacitance and ionic density profiles of a Coulomb fluid, both in the electrolyte and ionic liquid regimes. The calculation of electrostatic energy was performed using periodic Green’s functions. This approach allowed us to easily calculate the induced surface charge on the electrodes. The method requires only a small number of replicas to achieve any desired precision. To speed up the calculations, we precalculated the interaction potential at the beginning of simulations. The lattice gas formalism provides us with valuable insights, while maintaining simulation time feasible. We find that the transition between camel-shaped and bell-shaped differential capacitance curves depends on both the Bjerrum length and on the volume fraction of ionic liquid. For charge asymmetric systems, we find a secondary peak in the differential capacitance for large compacities. This peak is related to the appearance of additional structure in the double layer. In the strong coupling regime, the charge on the electrodes is overscreened—the first layer has charge which overcompensates the charge on the electrodes, resulting in a charge reversal. Finally, we find that the capacitance of size asymmetric ionic liquids is significantly smaller than that of symmetric ionic liquids, if all other parameters are the same. The great advantage of the lattice approach is that it significantly speeds up the simulations by allowing us to precalculate all the interactions. In the limit a → 0, at a fixed ionic size, we should recover the continuum limit. The crossover between lattice and continuum simulations will be the subject of future work.⁷¹ 

This work was partially supported by the CNPq, National Institute of Science and Technology Complex Fluids (INCT-FCx), CAPES, Alexander von Humboldt Foundation, and by the US-AFOSR under the Grant No. FA9550-16-1-0280.