Charge Relaxation Dynamics of an Electrolytic Nanocapacitor.

Understanding ion relaxation dynamics in overlapping electric double layers (EDLs) is critical for the development of efficient nanotechnology-based electrochemical energy storage, electrochemomechanical energy conversion, and bioelectrochemical sensing devices as well as the controlled synthesis of nanostructured materials. Here, a lattice Boltzmann (LB) method is employed to simulate an electrolytic nanocapacitor subjected to a step potential at t = 0 for various degrees of EDL overlap, solvent viscosities, ratios of cation-to-anion diffusivity, and electrode separations. The use of a novel continuously varying and Galilean-invariant molecular-speed-dependent relaxation time (MSDRT) with the LB equation recovers a correct microscopic description of the molecular-collision phenomena and enhances the stability of the LB algorithm. Results for large EDL overlaps indicated oscillatory behavior for the ionic current density, in contrast to monotonic relaxation to equilibrium for low EDL overlaps. Further, at low solvent viscosities and large EDL overlaps, anomalous plasmalike spatial oscillations of the electric field were observed that appeared to be purely an effect of nanoscale confinement. Employing MSDRT in our simulations enabled modeling of the fundamental physics of the transient charge relaxation dynamics in electrochemical systems operating away from equilibrium wherein Nernst-Einstein relation is known to be violated.

Introduction

Electric double layers (EDLs) associated with biomolecules, polymers, charged surfaces, and electrochemical interfaces play an important role in a wide variety of phase transformation,[1−3] interface, and transport phenomena[4−7] under the application of electric potentials or currents that vary over time. Advancing an understanding of transient ion transport and relaxation dynamics in EDLs is therefore critical not only for the development of nanotechnology-based efficient electrochemical energy storage,[8−10] electrochemomechanical energy conversion,[11] bioelectrochemical sensing,[12,13] and molecular trapping devices[14,15] but also for the synthesis of novel nanostructured materials through the control of EDL-mediated phase transformations[2,3] and electrochemical reactions.[16,17] Over the years EDLs have been modeled using equivalent circuit models,[18,19] analytic solutions, simulations of the Poisson–Nernst–Planck system of equations,[20−23] and, more recently, ab initio molecular dynamics (MD) and Monte Carlo (MC) methods.[24−26] Equivalent circuit models for EDLs are clearly oversimplified and suffer from drawbacks,[22,27,28] whereas the analytical or numerical solutions of the Poisson–Nernst–Planck system of equations that describe the dynamics of EDLs assume continuum dynamics.[29,30] Thus, in mesoscale systems with overlapping EDLs and relatively large Knudsen numbers (ratios of the molecular mean free path to the characteristic system length scale) the use of the Poisson–Nernst–Planck approach becomes untenable.[29,30] MD and MC methods have been successful in predicting equilibrium EDL structure and nonequilibrium steady-state properties of aqueous electrolytes,[24−26] but, using MD and MC methods, the simulation of the transient dynamics of charge relaxation in overlapping EDLs resulting from applied electric potentials or currents has remained a challenge. This is primarily because of the sheer computational expense required to (a) track individual ions and solvent molecules and (b) achieve statistical accuracy for electrolytic solutions where the ionic number densities are several orders of magnitude smaller than those for the solvent.[25,31] In this context, the use of a computationally efficient mesoscale simulation technique, such as the lattice Boltzmann method (LBM) that is based on the Boltzmann transport equation and is valid at relatively large Knudsen numbers, holds promise.[30]

Lattice Boltzmann Equation and the Nernst–Einstein Relation

In LBM, the time evolution of the single-particle distribution function f of the species of the particles of type i in an electrolyte mixture can be described using the nondimensional discretized lattice Boltzmann equation (LBE),[32−34] represented aswhere t is the time, Δt is the time step, r is the space coordinate, τ is the relaxation time, eα is the discrete particle velocity, feq is the local equilibrium Maxwellian distribution, F is the discretized external force term, and i (s, Cn, or An) corresponds to solvent particles, cations, or anions, respectively.

So far, however, application of LBM based modeling to transient simulations of EDLs in electrochemical systems and electroosmotic flows has not been attempted and assumptions of bulk electroneutrality and thin EDLs, or fully developed ionic concentration and electrostatic potential profiles for overlapping EDLs, as initial conditions have been employed to obtain nonequilibrium steady-state solutions.[4,35−38] This is true in part because all of these studies employ single relaxation times based on constant bulk ionic diffusion coefficients in the LBE for the simulation of ionic electrodiffusion. The use of a single relaxation time in the LBE in the Bhatnagar–Gross–Krook (BGK) approximation,[39] corresponding to diffusion in the bulk electrolyte, implicitly assumes the validity of the Nernst–Einstein relation between the equilibrium ionic diffusion coefficients Deq and their mobilities μeq. However, because the Nernst–Einstein relation (Deq = μeqkT/ze) is derived purely from equilibrium considerations, it is well established through both experiments and simulations that the Nernst–Einstein relation is not valid for the dynamics of charged particles in systems away from equilibrium.[40−42] As a result, any attempt to simulate the transient dynamics of such highly nonlinear systems in the BGK approximation results in numerical instabilities. Thus, to account for deviations from the Nernst–Einstein relation that affect relaxation times for ionic distribution functions in the LBE and to simulate the transient dynamics of EDLs in electrochemical systems away from equilibrium, we propose here a molecular-speed-dependent relaxation time (MSDRT) based on Tait’s theory of the mean free path.[43,44] The proposed relaxation time, which is also more generally applicable to hydrodynamic flow simulations, is first shown to recover a correct microscopic description of collision phenomena for the classic entrance flow problem of plane Poiseuille flow in a channel.[45] The MSDRT approach is then employed to simulate and to understand the nature of transient charge-relaxation dynamics in the overlapping EDLs of an electrolytic nanocapacitor.

Methods

For all simulation results presented herein, a 2D nine-velocity (D2Q9) model of LBM was employed. The simulations were carried out using message-passing interface (MPICH 2.0)-based parallel codes implemented in FORTRAN; postprocessing and visualization of the results were carried out using MATLAB (The MathWorks, Natick, MA). The discretized equilibrium distribution functions used in the LBE were computed using the following equation:with α → 1–9; weights wα (w1 = 4/9, w2–5 = 1/9, and w6–9 = 1/36); discrete particle velocities eα (e1 = [0, 0], e2,4 = [ ±c, 0], e3,5 = [0, ±c], e6,7 = [ ±c, c], and e8,9 = [ ±c, −c]); the lattice speed of sound cs = Δx/Δt = √3cs; and the speed of sound in the fluid medium represented as cs = (RT)1/2 = 1/√3. It is worth noting here that the LBE recovers the advection–diffusion equation for the ions in the hydrodynamic limit.[37]

Molecular-Speed-Dependent Relaxation Time

Typically, in the BGK approximation, the dimensionless relaxation time in the LBE is related to the equilibrium bulk kinematic viscosity of the solvent νseq or the ionic diffusion coefficients Deq, represented aswhere [νs, D → νseq, Deq]; Δx represents the lattice spacing; R is the gas constant; and T is the temperature. On the basis of kinetic theory,[46] in terms of the equilibrium mean free path leq for the particles in the electrolyte mixture, the equilibrium ionic diffusion coefficients and solvent viscosity can be expressed as Deq, νseq = (1/3)c̅leq, where c̅ = (8kT/πm)1/2 is the mean thermal speed of the particles of molecular mass m at equilibrium, T is the system temperature, and k is the Boltzmann constant. Away from equilibrium, the net macroscopic velocity of the particles u(r, t) also represents a statistical average of the microscopic velocity of the fluid particles. As such, at any given time and position in the simulation domain, the average local speed c(r, t) of the particles can be expressed as a sum of the molecular speed c̅ caused by random Brownian motion and the magnitude of the macroscopic velocity u(r, t), i.e.,Thus, to compute the MSDRT for systems away from equilibrium, we define the local kinematic viscosity νs(r, t) or the ionic diffusion coefficient D(r, t) aswhere l(c) is the local nonequilibrium mean free path. In the LB algorithm, because the local macroscopic velocity u(r, t) (resulting from the action of the external body force F on the i-th type of particle) is easily computed using equations for the macroscopic particle density n and the momentum density represented bythe problem of computing an MSDRT using eqs 3–7 is reduced to obtaining an estimate of the local nonequilibrium mean free path l that depends on the local speed c of the particles. This can be accomplished on the basis of Tait’s theory[44] with the underlying assumption that there are sufficient numbers of collisions in the system to allow for the existence of a local Maxwellian distribution of particle speeds. The number of collisions occurring in an interval of time dt between pairs of molecules of types i and j (with corresponding masses m and m in a mixture) can be expressed as[43]such that c, c, ψ, and ϵ lie in the ranges dc, dc, dψ, and dϵ, respectively, in the neighborhoods of c, c, ψ, and ϵ, respectively. Herein, c and c are the particle velocities, cr is the magnitude of the relative velocity (c = c – c) of the particles, ψ (= (π – χ)/2) is the angle between the relative velocity c and the unit vector k joining the centers of mass of particles i and j at the point of their closest approach, χ is the angle of deflection between the relative velocity c (before collision) and c′ (after collision), d is the average diameter of the molecules of types i and j, ϵ is the angle between plane LMN and the plane containing AP′ and the fixed z axis O (for a description of the geometry of collision, see Figure 1), and ξ and ξ are the Maxwell–Boltzmann molecular velocity distribution functions represented aswhere η is the number density of the particles of type i. Now, the total number of collisions during dt such that c lies in the velocity-space volume element d centered at is proportional to dt and to the number of particles ξ dc of type i, which can be expressed aswhere γ(c) is the collision frequency for molecules of type i moving at speed c with molecules of type j. The collision frequency γ(c) is independent of the direction of c and can be obtained by dividing eq 8 by ξ dc dt and integrating it over all values of c, ψ, and ϵ, represented as[43]where p = c(m/2kT)1/2, erf denotes the unnormalized error function, k is the Boltzmann constant, and d is the average molecular diameter of the molecules of type i and j. Using eq 11, we calculated the mean free path of particles of type i moving at speed c in the fluid mixture as l(c) = c/(Σ γ), which can then be substituted into eqs 5 and 3 to calculate the local nonequilibrium MSDRTs for use with the LBE. Because γ(c) and l(c) are both independent of the direction of c, the local nonequilibrium MSDRTs thus calculated are isotropic.

Figure 1

Geometry of a binary collision for particles of types i (green circle) and j (orange circle).

On the basis of the dependence of the MSDRT on the sum of the local mean thermal speed and the macroscopic speed of the particles, it may mistakenly be inferred that this would result in a violation of Galilean invariance. However, that is not true because the mean thermal speed depends only on the distribution of particle velocities irrespective of the inertial frame of reference employed and is therefore Galilean invariant[43] while the local macroscopic particle speed complies with Galilean invariance within an error on the order of the square of the Mach number O(Ma2) for the D2Q9 model employed here.[47] Therefore, the proposed LB formalism employing MSDRT is Galilean-invariant for the small particle velocities observed in our simulation results.

Entrance Flow Problem

Free-slip boundary conditions were specified from the start of the simulation domain to the inlet of the channel for a distance of 100 lu (lattice units) to simulate the effect of the reservoir, whereas no-slip boundary conditions were specified on the channel walls for the rest of the simulation domain length, i.e., 700 lu. The free-slip boundary conditions were implemented using specular reflection of the density distribution functions fs,α, such that there was no friction exerted on the fluid flow and the tangential component of the momentum remained unchanged.[48] The no-slip boundary conditions, designed to implement zero tangential fluid flux along the channel walls, were implemented using bounce-back of the particle distribution functions.[48] Using nonequilibrium bounce-back of the density distribution functions normal to the inlet and outlet boundaries,[49] a velocity boundary condition of us = [0.05,0] lu was imposed at the free-slip inlet while at the outlet a constant pressure boundary condition was imposed by fixing the fluid density at ns,out = 1 lu. The dimensional mass of the fluid particles was specified as ms = 18.015 amu. Other variables were set as follows: temperature as T = 298 K, density as ns = 1000 kg/m3, and diameter as ds = 1.92 Å. The theoretical fully developed Poiseuille flow velocity profile was computed using u = 4umax(y/h)(1 – y/h), where h is the distance between the channel walls, whereas the nondimensional pressure was computed from the local fluid density using the equation of state for the fluid p = ρRT.

Electrolytic Nanocapacitor

The time evolution of the ionic density distribution functions in response to the applied electric potential on the nanocapacitor electrodes was simulated using the LBE (eq 1) coupled to the nonlinear Poisson–Boltzmann equation. In our simulations, the electric field and potential, corresponding to the evolving ionic density distribution functions solved for using the LBE at each time step, were computed using the lattice Poisson–Boltzmann method (LPBM) proposed by Wang et al.[36] for the solution of the nonlinear Poisson–Boltzmann equation. Our choice of LPBM for the solution of the nonlinear Poisson–Boltzmann equation was influenced by its ease of parallelization. However, other more efficient multigrid Poisson–Boltzmann solvers[50,51] may also be employed equivalently. Dirichlet boundary conditions were applied for the electric potential Ve on the electrode surfaces such that Ve(x, y = 0) = +10 mV and Ve(x, y = h) = −10 mV for t ≥ 0, whereas an initial condition of a linearly varying electric potential (from +10 mV at the cathode to −10 mV at the anode) was specified throughout the simulation domain. To simulate perfectly blocking electrodes, we applied no-flux boundary conditions for the ions on the electrode surfaces using eq 7. On the simulation domain boundaries perpendicular to the electrode surfaces, we implemented periodic boundary conditions for the distribution functions used in the LPBM for the computation of the electric potential and in the LBE for the simulation of ionic drift–diffusion. To account for the electric and viscous drag forces acting on the ions in the solvent, we specified the discrete external force term F in eq 1 as[33,34]where E is the electric field, μs is the dynamic viscosity of the solvent, and ϕ is the number fraction of the ions of type i in the electrolyte mixture. The first and the second terms in the square brackets in eq 12 correspond to the electric and the viscous drag forces acting on the ions, respectively. For electrode spacings of h = 50 and 100 nm at t = 0, ions in the electrolyte experienced initial electric fields of E = 4 × 105 and 2 × 105 V/m, respectively. To understand the relaxation dynamics of the ions in the region between the two electrodes, simulations were carried out for various degrees of EDL overlap α (≡ 10 λD/h, where λD is the Debye length given by λD = (εkT/2nbz2e2)1/2 and nb is the bulk ionic concentration), solvent viscosity μs, electrode separations of 50 and 100 nm, and cation-to-anion diffusion coefficient ratios of 1:1 and 2:1. The equilibrium anion diffusion coefficient was specified as DAneq = 4 × 10–9 or 2 × 10–9 m2/s, depending on the diffusion coefficient ratio, whereas the equilibrium cation diffusion coefficient was fixed at DCneq = 4 × 10–9 m2/s. Because the electrode spacing for almost all simulations was fixed at h = 50 nm, the EDL overlap parameter was changed to α = 1.4, 0.9, 0.65, 0.2, or 0.1 by varying the bulk electrolyte concentrations (nb = 1.89, 4.57, 8.76, 92.54, or 370.19 mM, respectively). The effect of electrode spacing was studied by changing h to 100 nm and keeping the electrolyte concentrations fixed. The temperature for all simulations was fixed at T = 298 K; the solvent molecular mass ms and density ns were fixed at 18.015 amu and 1000 kg/m3, respectively; solvent relative dielectric constant was set to 78.547; and molecular masses of hydrated ions were fixed at 113.065 amu, corresponding to 5 molecular masses of the solvent molecule in the hydration shell and an ionic mass of mion = 22.99 amu. Corresponding to these conditions, the ionic diameters for the two anion to cation diffusion coefficient ratios of 1:1 and 2:1 were obtained from the Newton–Raphson method based iterative solution of the set of coupled quadratic equations (represented by eq 11) as dCn = dAn ≈ 3.48 Å and dCn ≈ 3.48 Å, dAn ≈ 5.71 Å, respectively. A lattice spacing of Δx = 0.5 Å corresponding to a simulation time step of Δt = 0.58 ps was employed for all simulations.

Results and Discussion

The entrance flow problem in isothermal hydrodynamics relates to the development of a fluid-flow profile in a channel that connects two reservoirs with unequal hydrostatic pressures.[45] It is assumed that as the fluid flows between the two reservoirs over time no change is observed in the pressure at either the inlet or the outlet of the channel. In the steady state, a parabolic flow velocity profile, corresponding to Poiseuille flow, develops in the channel. To validate the proposed MSDRT, we simulated the dynamics of a pure fluid flow in the absence of an external body force field using eq 1 until steady state was reached.

The simulation results in Figure 2a indicate the steady-state pressure along the length of the channel. The hydrodynamic pressure close to the channel walls at the point of transition (100 lu from a free-slip to a no-slip boundary condition for the fluid flow) is higher than it is in the center of the channel at the same position perpendicular to the channel wall. But despite the pressure being higher in a region close to the point of transition in boundary conditions, the fluid close to the walls would be forced to relax at the same rate as the fluid in the center of the channel if the BGK approximation is employed. A physically correct description of microscopic collision phenomena based on kinetic theory, however, must have a higher rate of relaxation for the fluid in a region of higher pressure on account of greater density than in a region of lower pressure. A survey of Figure 2a–c clearly shows that as expected with the use of the MSDRT in the LBE the nondimensional relaxation frequency ωs(r, t) is elevated in regions of higher pressure and lower speed and vice versa. Finally, Figure 2d illustrates how the transverse flow speed profile across the width of the channel at the outlet exactly matches that of the theoretically obtained parabolic profile for the plane Poiseuille flow.

Figure 2

Entrance flow in a channel. Nondimensional data in lattice units (lu): (a) pressure, (b) flow speed, (c) collision frequency, and (d) Poiseuille flow velocity profile across the channel width at the outlet of the channel. In a, the vertical arrow marks the point of transition from the free-slip boundary condition on the channel walls (on the left) to the no-slip boundary condition on the channel walls (on the right).

Electrolytic Nanocapacitor

For the simulation of the charging dynamics of the electrolytic nanocapacitor, we considered a primitive model of a symmetric 1:1 electrolyte. The ions in the solvent were described as point particles, and the effect of the solvent was included in the simulations by describing it as a background medium with a characteristic dielectric constant, constant density, and viscosity. Although the LB algorithm presented here can be easily extended to simulate coupled ion–solvent relaxation dynamics as well, a deliberate choice was made to exclude the dynamics of the solvent from the simulations for the sake of simplicity and therefore the solvent velocity us was approximated to zero. However, for the purpose of computing the MSDRTs for the ionic distribution functions f, hydrated ions with an effective mass m and finite diameter d were considered for use in eq 11. The time evolution of ionic density distributions in the overlapping EDLs of the electrolytic nanocapacitor in response to an applied electric potential step at t = 0 was simulated using the discrete LBE (eq 1). The resulting electric force, acting on the uniformly distributed electrolyte ions at t = 0, caused the ions to accelerate towards electrodes of opposite polarities and impinge on the capacitor electrodes. The subsequent ion relaxation dynamics in the region between the two electrodes were then investigated.

Space-Averaged Current Density

The effect of EDL overlap α and solvent viscosity on the space-averaged current densityis depicted in Figure 3a–c. It is observed that as the degree of EDL overlap α is reduced from 1.4 to 0.65 and then to 0.1, the behavior of the space-averaged current density changes gradually from oscillatory to monotonic accompanied by a decrease in the maximum amplitude. Increasing the solvent viscosity μs from 0.0008 to 0.0018 Pa·s dampens the oscillatory behavior for α = 1.4, 0.65 and reduces the amplitude of the oscillations (Figure 3a,b); in the case of α = 0.1, however, it increases the relaxation time to equilibrium for the space-averaged current density (Figure 3c). The oscillatory behavior in the case of α = 1.4, 0.65 is seen to persist for μs = 0.0008, 0.0013 Pa·s much longer than what might be expected (Figure 2a,b). Figure 3d,e shows that increasing the electrode separation h from 50 to 100 nm decreased the maximum amplitude of Javg for both the bulk ionic concentrations of nb = 4.6 and 92.5 mM. Doubling the electrode separation also increased the period of oscillations in Javg by a factor of ∼1.5–2.0 for nb = 4.6 mM and doubled the relaxation time to equilibrium for nb = 92.5 mM. The effect of the ratio of equilibrium cation-to-anion diffusion coefficients on Javg for a solvent viscosity of 0.0018 Pa·s is shown in Figure 3f,g. Reducing the anion diffusivity by half reduced the peak Javg for both α = 1.4 and 0.1 and was accompanied by a slight increase in the relaxation time for α = 0.1 (Figure 3f,g). The disparity in the diffusivity of ions also introduced persistent oscillations with very small amplitudes in Javg for α = 1.4 versus those observed when the ionic diffusivities were equal.

Figure 3

Space-averaged current density as a function of time. Effect of EDL overlap and solvent viscosity with h = 50 nm, DCn/DAn = 1:1, and α = (a) 1.40, (b) 0.65, and (c) 0.10; electrode separation with μs = 0.0013 Pa·s, DCn/DAn = 1:1, and nb = (d) 4.6 mM and (e) 92.5 mM; ratio of cation-to-anion diffusion coefficient with μs = 0.0018 Pa·s, h = 50 nm, and α = (f) 1.4 and (g) 0.10.

In terms of the underlying physics, the decrease in the peak amplitude of Javg with an increase in solvent viscosity μs (Figure 3a–c) and a halving of the anion diffusion coefficient DAn (Figure 3f,g) are similar and result from a reduction in the mobility of ions. This is true because both viscosity and diffusivity are empirically related to each other and mobility through the Nernst–Einstein and Stokes–Einstein relations D = μkT/ze = kT/3πdμs. A decrease in ionic mobility therefore leads to decreased conductivity of the electrolytic medium between the capacitor electrodes, resulting in an increase in the potential drop across the electrolyte and a consequent decrease in the peak amplitude of Javg. This also explains faster relaxation to equilibrium for the monotonic case of α = 0.1 with decreasing viscosity of the solvent as more mobile ions charge up the electrodes faster (Figure 3c). For α = 1.4 and 0.65, a decrease in viscosity and the consequent increase in ionic mobility cause the ions to oscillate longer (Figure 3a,b).

Now, for a monotonically charging capacitor, the current density as a function of time is J = J0e–, where τc = RC = (h/σA)C is the characteristic relaxation time for the capacitor, A is the effective area of the capacitor electrodes, C is the capacitance, σ is the conductivity of the dielectric medium between the two electrodes, and J0 is the peak current amplitude. Thus, a doubling of the electrode spacing h must also result in a doubling of the characteristic charging time τc (Figure 3e). An increase in the resistance R = (h/σA) also accounts for the decrease in the maximum value for Javg (Figure 3d,e). Similarly, doubling the distance between the capacitor electrodes doubles the time taken by ions to move from one electrode to another, thereby resulting in a near doubling of the period of oscillations for Javg (Figure 3d). For thin EDL overlaps, the results for Javg indicate that the monotonic charging dynamics for an electrolytic nanocapacitor is recovered (as in macroscopic electrochemical systems), whereas for large EDL overlaps, oscillatory behavior is observed.

Collective Oscillatory Dynamics of Ions

For an EDL overlap of α = 1.4, Figure 4a,b depicts the damped harmonic motion of ions as the ionic velocities periodically reverse direction between the perfectly blocking electrodes and then decay to zero over time. As a consequence, the electric potential and field profiles also evolve in an oscillatory fashion over time trending towards their respective equilibrium electrostatic profiles (Figure 4c,d and Movie S1). In the case of α = 0.1, Figure 5a–d and Movie S2 show a monotonic convergence to equilibrium electrostatic values for these variables. The force that ions experience at any point between the two electrodes of the capacitor is proportional to electric field E experienced by the ions at that location. A comparison of the electric field profiles in the gap between the capacitor electrodes in Figures 3d and 4d for the two cases of EDL overlap shows that both cases exhibit a trough in the middle. This trough, in the case of an EDL overlap of α = 1.4, is almost parabolic in shape with a nonzero electric field in the center, whereas for α = 0.1, the trough is flat-bottomed with a field that vanishes in the center of the capacitor. The oscillatory and monotonic motions of ions can be likened to the motion of a ball rolled in parabolic and flat-bottomed troughs, respectively, with frictional surfaces under the influence of gravity from the top edges of the troughs. In the parabolic case, the ball executes a back-and-forth oscillatory motion before coming to rest, whereas in the flat-bottomed case, it comes to rest more or less monotonically. In the case of overlapping EDLs, this indicates that the parabolic or near-parabolic shape of the electric field profile in the center of the nanocapacitor in conjunction with a nonvanishing electric field is responsible for the observed oscillatory behavior.

Figure 4

Oscillatory behavior of ions for large EDL overlap as a function of time for α = 1.4, μs = 0.0018 Pa·s, h = 50 nm, and DCn/DAn = 1:1 (see also Movie S1): (a) cation velocities, (b) anion velocities, (c) electric potential, and (d) electric field.

Figure 5

Monotonic charge relaxation for small EDL overlap as a function of time for α = 0.10, μs = 0.0018 Pa·s, h = 50 nm, and DCn/DAn = 1:1 (see also Movie S2): (a) cation velocities, (b) anion velocities, (c) electric potential, and (d) electric field.

Anomalous Spatial Oscillations of the Electric Field

A persistence of oscillations in Javg was observed for three cases with an EDL overlap of α = 1.4: (i) μs = 0.0008 Pa·s, (ii) μs = 0.0013 Pa·s (Figure 3a), and (iii) μs = 0.0018 Pa·s (when an asymmetry was introduced by changing the cation-to-anion diffusion coefficient ratio from 1:1 to 2:1, Figure 3f). Since the origin of persistent oscillations is similar for cases (i) and (ii), we further analyzed cases (i) and (iii) only.

Considering case (i) first, the ionic velocities plotted in Figure 6a,b show a general oscillatory behavior for a period up to 11 ns which is reflected in the dimensionless collision frequencies for the ions, electric field, the spatial cross-correlation (Jcorr) of ionic flux densities (JCnf and JAnf) and space-averaged current densities (Javg, JCnavg, and JAnavg) (Figure 6a–e and Supporting Information Movie S3). Between 11 and 20 ns, two nodes spontaneously develop at roughly one-third of the electrode spacing away from either electrode; at these nodes, the ionic velocities approach zero and then switch signs (Figure 6a,b). The region of these two nodes between the nanocapacitor electrodes is marked by intense collisions, illustrated by the presence of two sharp peaks at roughly one-third of the electrode spacing away from either electrode (Movie S3). Once the persistent oscillations are fully developed between 58 and 93 ns, both cations and anions self-organize, resulting in their velocities forming a single node between the capacitor plates. At this node, the ionic velocities switch signs at any given instance in time (Figure 7a–d and Movie S3), and the node’s position in the center of the nanocapacitor is again marked by intense collisions between ions (Figure 7a,b). On either side of this node, there exist regions where the collision frequencies are much lower and the ionic velocities are high. In the regions of low collisions, the ionic velocities switch signs at different instances in time. Next, upon examination of the spatial cross-correlation Jcorr of the y component of the ionic flux densities JCnf(t) and JAnf(t) (Figure 6d) in the range of t = 0–11 ns when JCnavg and JAnavg are in-phase (Figure 6e), one observes that the spatial cross-correlation Jcorr stays negative and the electric field amplitude for E exhibits oscillatory behavior about the equilibrium electrostatic electric field in the center of the simulation domain (Figure 6c). The negative values of Jcorr in Figure 6d are consistent with the in-phase behavior of JCnavg and JAnavg because for a given electric field at a certain point oppositely charged ions must move in opposite directions, resulting in a negative cross-correlation of their respective fluxes. Between 11 and 20 ns, however, a positive correlation between ionic flux densities JCnf(t) and JAnf(t) begins to manifest itself (Movie S3 and Figure 6d). For the spatial cross-correlation Jcorr in the latter out-of-phase range of t = 58–93 ns for JCnavg and JAnavg, positive values are observed with a peak at a spatial lag (a separation of about 25.3 nm, Figure 7d). Electric field E also starts to exhibit spatial oscillations in this time range (Figure 7c and Movie S3), in contrast to the simple oscillations in amplitude in the center of the nanocapacitor observed when JCnavg and JAnavg are in phase (Figure 6c). It is seen that the minima of these spatial oscillations are again separated by approximately 25.3 nm.

Figure 6

Oscillatory charge relaxation before the onset of plasmalike collective behavior as a function of time for α = 1.40, μs = 0.0008 Pa·s, h = 50 nm, and DCn/DAn = 1:1 (see also Movie S3): (a) cation velocities, (b) anion velocities, (c) electric field, (d) spatial cross-correlation of cation and anion flux densities, and (e) space-averaged total, cation, and anion current densities (inset: transition of space-averaged cation and anion current densities from in-phase to out-of-phase behavior marking the onset of plasmalike collective behavior).

Figure 7

Anomalous oscillations in electric field after the onset of plasmalike collective behavior as a function of time for α = 1.40, μs = 0.0008 Pa·s, h = 50 nm, and DCn/DAn = 1:1 (see also Movie S3): (a) cation collision frequencies, (b) anion collision frequencies, (c) electric field, and (d) spatial cross-correlation of cation and anion flux densities. (e) Dielectric slab in the nonuniform electric field of a parallel plate capacitor.

Although much weaker than in case (i), a similar positive spatial cross-correlation of ionic flux densities is also observed under the conditions in case (iii) (EDL overlap of α = 1.4, cation-to-anion diffusion coefficient ratio of 2:1, solvent viscosity of μs = 0.0018 Pa·s; see Movie S4). Since such a behavior is absent in the case of a 1:1 diffusion coefficient ratio under similar conditions, the effect of introduction of an asymmetry in the diffusion coefficients or mobility of the ions is to cause persistent spatial oscillations in electric field E and a positive correlation in the ionic flux densities JCnf and JAnf. In contrast to the three cases discussed above, in the case of an EDL overlap of α = 0.1, a cation-to-anion diffusion coefficient ratio of 1:1, and a solvent viscosity of μs = 0.0018 Pa·s (Movie S2), the linearly negative Jcorr monotonically approaches zero as the system marches to equilibrium, as illustrated by the monotonic time evolution of electric field E.

Thus, it becomes apparent that when JCnavg and JAnavg are out of phase the oppositely charged ions in the two opposite halves of the simulation domain exhibit positively correlated motions and form a sort of dynamic standing-wave pattern (Movie S3). It is this collective behavior that results in plasmalike spatial oscillations in the electric field. But then what is the reason behind such anomalous collective behavior? The reason becomes clear when one analyzes the behavior of a dielectric slab upon insertion into a capacitor (Figure 7e). A dielectric slab, when slowly introduced into a capacitor, experiences a non-uniform electric field that exerts a force = ∇(·) on it which tends to drive it into a region of higher electric field. Because the field inside a dielectric slab points in a direction opposite to that of the electric field between the capacitor plates, the insertion of a dielectric slab has the effect of lowering the net electric field in the region that it occupies. Now, consider that this dielectric slab is made up of oppositely charged ions interacting with each other to form temporary dipoles. With this, if the binding energy of the dipoles in the dielectric slab was less than the interaction energy of the individual ions (forming the dipole) with the external electric field of the capacitor then the dipoles will be pulled apart into individual ions and forced to move away (in the same direction) from the region of higher electric field as both types of ions are in a position that is far away from the equilibrium for the individual ions.

Consider next a set of conditions with a region of lower electric field and low ionic velocities accompanied by intense collisions in the center of the nanocapacitor that allows them to interact again to form dipoles. Now, if the field between the capacitor plates be non-uniform (as it is in the case of the electrolytic nanocapacitor), the dipolar interaction between the ions would again force them to move to a higher electric field (this time in the opposite direction) because of a residual momentum from their previous interaction with a region of high electric field. The net electric field in the region occupied by the ions interacting like dipoles would also be lowered until the time the ions again encounter a region of high electric field that is sufficiently strong to break them apart and cause a repeat of the cycle described above. Such repetitive behavior would thus result in the oscillatory behavior of the electric field in space as a consequence of the spatially correlated motion of oppositely charged ions similar to that observed in Figure 7a–d. Thus, it can be concluded that the persistent spatial oscillations of E, and JCnavg and JAnavg are most likely due to a dipolar interaction occurring between the ions that results in a plasmalike collective behavior.

The plasma frequency ωp for an electrolyte in terms of the Debye length λD is given as ωp = λD–1(kT/m*)1/2, where m* is the harmonic mean of the effective cation and anion masses.[39] An EDL overlap of α = 1.4, corresponding to an electrolyte concentration of nb = 1.89 mM, gives a plasma frequency of ωp ≅ 20.6 GHz, whereas the collision frequency for the ions in a solvent of density ns = 1000 kg/m3 is on the order of νc ≈ 5 THz. For such an electrolyte where the collision frequency is νc ≫ ωp, the current density is always in phase with the applied electric field, and a plasma oscillation cannot be excited.[52] However, even with an oscillation frequency of 1.85 GHz (ωp/10) for the in-phase ionic current densities JCnavg and JAnavg just before the onset of the out-of-phase behavior, plasmalike spatial oscillations are observed (Figure 6e inset). Such behavior is clearly anomalous; therefore, the results pertaining to plasmalike collective behavior presented here for overlapping EDLs in an electrolytic nanocapacitor appear to be purely an effect of nanoscale confinement.

Conclusions

The use of a single relaxation time under BGK approximation in the LBE greatly simplifies computation. However, in general, the use of a single relaxation time in the BGK approximation in LBM-based simulations is known to suffer from problems of numerical instability and spurious artifacts, especially in systems away from equilibrium that involve strong spatial variations in body forces, particle densities, or temperature.[4,53−57] Various attempts to solve these problems have been made using multiple relaxation time (MRT),[58] two relaxation time (TRT),[59] and entropic[60] LB models. These approaches employ either constant relaxation times based on bulk-transport coefficients or a variable over-relaxation parameter derived from them using entropic constraints. Because relaxation times in MRT and TRT LB models are based on constant bulk-transport coefficients, neither of these two methods (aimed at improving the numerical stability of LBM simulations) recovers a correct microscopic description of the collision phenomena. It is also unclear if the computationally expensive entropic LB models achieve this either. Other variants of the LB algorithm have explored the heuristic introduction of density dependence into the relaxation time (through a simple division of the constant BGK relaxation time by the position-dependent local fluid density)[61] and constant relaxation time approaches (based on tunable diffusivities for binary mixtures) through a separate consideration of the mutual and cross-collisions that is aimed at recovering the prescribed transport properties of a binary mixture.[62−64] The MSDRT proposed herein, however, allows for continuously varying relaxation times in a simulation domain that recovers a correct microscopic description of the collision phenomena while simultaneously retaining the simplicity of the LB algorithm in the BGK approximation. Because MSDRT incorporates the local macroscopic speed of the particles, their number densities, their molecular masses, and temperature in the computation of the relaxation time for the LBE, it may be better suited for the simulation of the dynamics of complex systems, such as interfaces in multiphase fluids[53,54] or mesoscale electrochemical or electrokinetic systems[6,7] with overlapping EDLs[4,55] that suffer from problems of numerical instabilities or spurious artifacts and involve strong spatial variations in body forces,[4,55] particle densities,[53,54,57] or temperature.[56] The effectiveness of the MSDRT proposed here is demonstrated through the transient simulations for the highly nonlinear charge-relaxation dynamics of an electrolytic nanocapacitor with spatially nonuniform force fields generated by inhomogeneous charge distributions in overlapping EDLs.

Additionally, we have successfully shown that the proposed LBM that uses the continuously varying MSDRT can be employed to simulate and to explore the fundamental physics of nonequilibrium charge-relaxation dynamics in mesoscale electrochemical systems wherein the Nernst–Einstein relation is known to be violated. So far, this had been one of the most challenging unsolved problems in nonequilibrium electrokinetics. Furthermore, our results predict the presence of hitherto unobserved phenomena of spatial oscillations in an electric field arising from the spontaneous dipolelike interaction of ions in response to a step potential applied across an electrolytic nanocapacitor under the conditions of large EDL overlap and low solvent viscosity.

For these reasons, we believe that the LBM presented here will open avenues for advancing an understanding of phase transformation, transport, and interface phenomena in electrochemical, electrokinetic, and micro or nanofluidic systems critical to the development of technologically advanced mesoscale devices and the synthesis of novel nanomaterials.