Charging dynamics of electrical double layers inside a cylindrical pore: predicting the effects of arbitrary pore size

Abstract

Porous electrodes are found in energy storage devices such as supercapacitors and pseudocapacitors. However, the effect of electrode-pore-size distribution on their energy storage properties remains unclear. Here, we develop a model for the charging of electrical double layers inside a cylindrical pore for arbitrary pore size. We assume small applied potentials and perform a regular perturbation analysis to predict the evolution of electrical potential and ion concentrations in both the radial and axial directions. We validate our perturbation model with direct numerical simulations of the Poisson–Nernst–Planck equations, and obtain quantitative agreement between the two approaches for small and moderate potentials. Our analysis yields two main characteristic features of arbitrary pore size: (i) a monotonic decrease of the charging timescale with an increase in relative pore size (pore size relative to Debye length); (ii) large potential changes for overlapping double layers in a thin transition region, which we approximate mathematically by a jump discontinuity. We quantify the contributions of electromigration and charge diffusion fluxes, which provide mechanistic insights into the dependence of charging timescale and capacitance on pore size. We develop a modified transmission circuit model that captures the effect of arbitrary pore size and demonstrate that a time-dependent transition-region resistor needs to be included in the circuit. We also derive phenomenological expressions for average effective capacitance and charging timescale as a function of pore-size distribution. We show that the capacitance and charging timescale increase with smaller average pore sizes and with smaller polydispersity, resulting in a gain of energy density at a constant power density. Overall, our results advance the mechanistic understanding of electrical-double-layer charging.

---

1 Introduction

Batteries and fuel cells are traditional porous-material-based electrochemical devices. Over the past 60 years, new devices such as supercapacitors¹ have started to emerge. They are comprised of porous electrodes – typically made up of dispersions of activated carbon spheres – immersed in aqueous, organic, or ionic liquid electrolytes. The electrodes store charge through the physical adsorption of dissociated ions onto their pore surfaces, forming a charged region commonly referred to as the electrical double layer.^2,3 Obviating the need for redox reactions, these devices charge faster than batteries and present better cyclability. Nevertheless, their energy density and capacitance are limited by their available specific surface area.³ In view of these characteristics, supercapacitors bridge the gap between traditional capacitors and batteries, being used in situations where fast response and moderate energy output are required, e.g., stabilization of energy fluctuations in power grids⁴ and memory protection in electronic devices.⁵ More recently, hybrid capacitors have been designed in an effort to utilize the energy storage mechanisms of both electrical double layers and reduction-oxidation (redox) reactions.^6,7 They consist of two distinct electrodes, one supercapacitor-like which stores charge physically into double layers, and another metal oxide pseudocapacitor-like which accumulates energy by performing fast oxidation surface reactions.⁶ While there have been significant advances in the material design of supercapacitors and hybrid capacitor electrodes, the effect of pore-size distribution on the energy storage properties of these devices remains unclear.⁶

The earliest models for electrode charging address the geometry of an electrolyte between flat plates, dating back to the works of Gouy,⁸ with later contributions from Chapman⁹ and Stern¹⁰ to consolidate the well known Gouy–Chapman–Stern model. Notably, a wealth of different effects have been studied to complete this simplified geometrical picture of flat plate charging. Bazant et al.¹¹ reviewed the previous flat plate charging studies and developed an approach to solve the Poisson–Nernst–Planck (PNP) equations, from the linear regime at low voltages to the nonlinear effects at high voltages. Feicht et al.¹² compared one-dimensional PNP solutions with the experimentally observed reverse peaks in electrolytic cell discharging. Kilic and Bazant^13,14 derived modified PNP equations including ion-radius effects to study the influence of ion crowding in charge storage. For higher electrolyte concentrations, other effects such as multi-component diffusion given by Stefan–Maxwell fluxes,¹⁵ and ion correlations^16–18 have been studied, for instance in some continuum models.^19–21

The porous geometry is incorporated in the form of equivalent circuit representations, widely used in the modeling of pore charging in supercapacitor electrodes,^22–26 stemming from the pioneering work of de Levie.^27,28 Later, Biesheuvel and Bazant² extended the circuit to high potentials for capacitive deionization applications. Recent papers^26,29 further discuss the relationship between the equivalent circuit and the corresponding transmission line continuous equation for pores with finite lengths. Throughout most pore charging models, common assumptions are either of thin double layers^{2,27,28,30,63} within the pores, i.e., such that the length of the charged region is much smaller than the pore size, or overlapping double layers,^31,32 where the charged regions extending from the opposite sides of the surface meet. However, pore sizes can range from less than 2 nm for micropores to more than 50 nm for macropores.^6,25 Thus, a first-principles approach that extends pore charging models to arbitrary pore sizes is required in order to accurately describe the charging of supercapacitors and predict their properties, such as capacitance, energy density, and power density. Some of the works which address pore-size dependence focus on the equilibrium response,³³ while others propose transient transport-equation-based numerical schemes that describe porous network charging for arbitrary pore sizes.^34,35 However, to the best of our knowledge, the relative importance of charge transport mechanisms for arbitrary pore sizes and their implications on transmission line circuits have not been reported.

In our previous work,³⁶ an analytical model based on the PNP equations was proposed to describe the charging of pores at low applied potentials in the limit of overlapping double layers. Here, inspired by perturbation models on electrokinetics,^37–42 we develop a perturbation expansion model for arbitrary pore sizes – i.e., pore radii – in the limit of small applied potentials. We compare the predictions of the model to direct numerical simulations to show that the perturbation solution yields good results even for moderate applied potentials (≈50 mV). We demonstrate that a modified transmission line circuit (compared to classic supercapacitor literature^27,28) includes a resistor representing finite changes in electric potential at a thin entrance region at the mouth of a pore. We derive an effective capacitance to study the effect of pore-size distribution on energy and power densities. Besides addressing transient charging, the solution developed here quantifies the contributions of diffusion and electromigration and provides mechanistic insights into the charging process of pores of arbitrary sizes.

2 Problem formulation

We consider an ideally conducting porous electrode that consists of tortuous pores of different radii, filled with an electrolyte; see the schematic in Fig. 1a. Once a voltage difference is applied, counterions are attracted to the pore surfaces, forming electrical double layers of thicknesses that may be thin or comparable to the pore size.^3,43 We follow common practice^2,11,36 in assuming the formation of a static diffusion layer (SDL), an electroneutral region beside the electrode. Fig. 1b shows a simplified setup with a cylindrical pore of length [small script l]_p and radius a_p, where the radial and axial directions are denoted by r and z, respectively. The SDL, of length [small script l]_s and radius a_s, is adjacent to the pore. We denote the cation and anion concentrations by c_±(r, z, t) and electric potential by ψ(r, z, t), where t is time. The position z = 0 represents the interface between the pore and the SDL. We assume that the SDL is in contact with the bulk such that c_±(r, −[small script l]_s, t) = c₀ and ψ(r, −[small script l]_s, t) = 0. At t = 0, the concentration of ions everywhere is c_±(r, z, 0) = c₀. At t > 0, since the electrode material is an ideal conductor, the potential at the surface of the pore ψ(a_p, z, t) = ψ_D. We also assume that the pore surface is ideally blocking, i.e., the flux of ions across the pore surface is zero. In addition, we assume that , , the ions are monovalent, that ion diffusivities inside the pore are equal and given by D_p. The ion diffusivities outside the pore are also assumed equal and are given by D_s. We note that the ion diffusivities inside and outside the pore may be different due to confinement effects. As we show later, the ratio of diffusivities – i.e., – dictates the interaction between the pore and the SDL. We denote Debye length by λ such that , where ε is the electrical permittivity. For reference, a 1 M aqueous electrolyte at room temperature has a Debye length λ ≈ 0.3 nm. The objective of this article is to determine the value of ψ(r, z, t) and c_±(r, z, t) for an arbitrary relative pore size .

Fig. 1 A schematic of the problem setup. (a) A porous electrode is subjected to an applied voltage. The electrode pores are filled with an electrolyte with cations (in red) and anions (in blue). The pore electrolyte is in contact with a static diffusion layer. (b) Schematic representation of a single cylindrical pore (using dimensionless variables) and the division of the domain into bulk, static diffusion layer (SDL), transition, and pore regions. The bulk has Ψ = 0 and equal concentration of both ionic species c_± = 1. The SDL radial boundary, with vanishing normal gradients, is represented by dashed lines. The potential at the surface of the pore is Ψ_D, the lengths of the static diffusion layer and pore are [small script l]_s/[small script l]_p and 1, respectively. The radii of SDL and the pore are a_s/a_p and 1, respectively, and the Debye length is λ/a_p.

Physically, when the potential is applied on the surface of a pore, oppositely charged ions migrate inside the pore, while the similarly charged ions transport out of the pore. This relative transport of ions produces an electrical current. In addition, due to the charge imbalance, an electrical double layer forms adjacent to the surface. As time progresses, the potential drop across the electrical double layer starts to saturate, ions stop migrating, and the current ceases.

To solve for ψ and c_±, we start by writing the Poisson–Nernst–Planck equations^44–46

−ε∇²ψ = e(c₊ − c₋), (1b)

where N_± are the ion fluxes. Inside the poreand outside the poreWe note that N_± = e_rN_±r + e_zN_±z by angular symmetry. We introduce dimensionless charge density , salt concentration , potential , time , axial position , radial position , and gradient operator . Note that . In addition, 0 ≤ R ≤ 1 for 0 ≤ Z ≤ 1 (i.e., the pore region), and for (i.e., the static diffusion layer region; see Fig. 1b). With these substitutions, the set of eqn (1) becomeswhere , is the dimensionless charge flux and is the dimensionless salt flux. Inside the poreand outside the poreAt τ = 0⁺, inside the pore, the potential is constant and equal to the wall potential since the electrical double layer hasn't developed yet. In contrast, the potential is linear in the SDL due to electroneutrality. Therefore, eqn (1) are subjected to the following initial conditions

ρ(R, Z, 0⁺) = 0, (3a)

s(R, Z, 0⁺) = 2, (3b)

Ψ(R, Z ≥ 0, 0⁺) = Ψ_D. (3d)

The boundary conditions (BCs) at are given by the bulk condition, or

At the end of the pore, i.e., Z = 1, the BCs are that gradients of potential and concentration vanish for l_p ≫ a_p, implying

Due to the symmetry, the BCs at the center of the system, i.e., R = 0, are simply

At the surface of the pore, i.e., R = 1 and Z ≥ 0, the BCs are of ideally blocking electrode and constant potential, which are given by

J_R|_R=1 = 0, (7a)

W_R|_R=1 = 0, (7b)

Ψ(1, Z > 0, τ) = Ψ_D. (7c)

We note that eqn (1c) represents the ideally conducting electrode BC. Finally, at the boundaries of the SDL, i.e., and Z < 0, the BCs are vanishing gradients, orThis assumption is valid for non-interacting pores where the radial currents are identically zero, consistent with the treatment of SDL in the literature.^2,29,36 We solve the set of eqn (2)–(8) numerically using OpenFOAM.^47,48 The details of geometry, mesh, and algorithm have been described in ref. 36. We refer to the solution from OpenFOAM as direct numerical simulations (DNS). Next, we perform a regular perturbation analysis⁴⁴ to obtain an analytical expression for ρ(R, Z, τ) and Ψ(R, Z, τ).

3 Regular perturbation analysis

In this section, we focus on the small potential limit, i.e., |Ψ_D| ≪ 1, and conduct a regular perturbation analysis to describe the charge and potential inside the pore. To this end, we divide the solution into three regions: (I) SDL, (II) inside the pore, and (III) transition region between the SDL and the mouth of the pore.

3.1 Static diffusion layer

The SDL is characterized by electroneutrality, i.e., ρ = 0. Furthermore, for |Ψ_D| ≪ 1, s = 2.^2,36 Therefore, as per eqn (2c) and the boundary conditions in eqn (4c), (6), and (8), it is straightforward to show that Ψ is linear in Z and is independent of R. Mathematically, we writewhere Ψ_left = Ψ(R, 0⁻, τ), i.e., the centerline potential to the left of the SDL–pore interface, and is to be determined.

3.2 Inside the pore

The region inside the pore consists of ion concentration and potential varying in time as well as in both radial and axial directions. In the low-potential limit, we propose regular perturbation expansions of the dependent variables in the small parameter Ψ_D, i.e.,

ρ = ρ₀ + ρ₁Ψ_D + O(Ψ_D²), (10a)

s = s₀ + s₁Ψ_D + O(Ψ_D²), (10b)

Ψ = Ψ₀ + Ψ₁Ψ_D + O(Ψ_D²). (10c)

We also introduce variables ρ_m(Z, τ) and Ψ_m(Z, τ) that represent the centerline charge and density. As per our proposed perturbation expansion,

ρ_m = ρ_m0 + ρ_m1Ψ_D + O(Ψ_D²) (10d)

Ψ_m = Ψ_m0 + Ψ_m1Ψ_D + O(Ψ_D²). (10e)

In this article, we only focus on the leading-order and first-order terms.

The leading-order terms, i.e., ρ₀, s₀ and Ψ₀ are obtained from the response in the absence of an applied potential, i.e., Ψ_D = 0. Thus, it can be seen that the leading-order coefficients ρ₀ = 0, s₀ = 2 and Ψ₀ = 0 satisfy the set of eqn (2) with the initial conditions (3) and boundary conditions (5)–(7).

Inserting the expansions provided in eqn (10) (with ρ₀ = 0, s = 2 and Ψ₀ = 0) in eqn (2a) and (2d) and collecting the first-order terms in Ψ_D yields

Similarly, Poisson's equation – eqn (1b) – after substitution of the perturbation expansions takes the formOnly the known leading-order term in the salt concentration enters eqn (11), such that eqn (11) and (12) suffice to determine first-order corrections to charge density and potential in the low potential regime. Therefore, we do not need to further solve eqn (2b).

Since , the diffusion in the radial direction is much faster than in the axial direction and we can assume quasiequilibrium in the radial direction.³⁴ Thus, with symmetry and ideally blocking conditions (see eqn (6) and (7)), we have J_R = 0. Utilizing this condition in eqn (2d) and collecting the first-order terms in the perturbation expansion, we get

which implies that³⁶

ρ₁ + 2Ψ₁ = ρ_m1 + 2Ψ_m1, (14)

where ρ_m1(Z, τ) and Ψ_m1(Z, τ) are to be determined. Next, by utilizing in eqn (12), we get³⁶Integrating eqn (15) with and Ψ₁(Z, 1, τ) = 1 yieldswhere I_n is the modified Bessel function of the first kind of order n. By substituting R = 0 in eqn (16), we obtainNext, by combining eqn (16) and (17), substituting the leading- and first-order coefficients in eqn (10), and neglecting higher-order terms, we getSimilarly, substituting eqn (16) and (17) into eqn (14) and neglecting higher-order terms in the perturbation expansion of ρ, we getNext, we note that since J_R = 0, so we can simplify eqn (11) to obtainTo determine the axial dependence, we average eqn (20) over the cross-sectional area of the pore (integrating across the radial direction) by utilizing the known radial dependence in eqn (18) and (19). Mathematically,We emphasize that this crucial step enables us to derive a solution for arbitrary , thereby bridging the previously reported trends in thin and overlapping double layer limits.^2,27,28,36 Now, substituting eqn (18) and (19) into eqn (21), we obtainThe boundary conditions for eqn (22) include Ψ_right = Ψ_m(0⁺, τ) and , where Ψ_right, i.e., the centerline potential to the right of the SDL-pore interface, is to be determined.

3.3 Transition region

We now focus on the transition region between the SDL and the mouth of the pore. First, we emphasize that the transition region is only relevant for overlapping double layer limits. In the thin double layer limit, ρ = 0 inside the majority of the pore as well as the SDL region, and thus the transition region becomes irrelevant. This is consistent with the derivations in Biesheuvel et al.,^2,30 who assume a continuous variation in Ψ across the SDL–pore interface in the thin double layer limit, thereby implicitly neglecting the transition region.

In contrast, in the overlapping double layer limit, ρ inside the pore is non-zero whereas ρ in the SDL region is zero, and a transition region is required to connect the two regions. Accordingly, inside the transition region, ρ varies from zero to the value inside the pore. To estimate the thickness of the transition region δ (scaled by [small script l]_p), we emphasize that charge density is related to the length scale of the charge gradients by the Poisson equation; see eqn (2c). Because the radial variation inside the SDL can be ignored, the dimensionless length scale over which charge gradients could be present inside the SDL is . Therefore, using eqn (2c) and assuming , it is straightforward to show that , enabling us to assume electroneutrality inside the SDL. In contrast, the smallest length scale over which the charge gradients are present inside the pore is , which implies ;^{36,44,49–51} see also eqn (29). Next, in the transition region, ρ varies from zero to the value inside the pore. Therefore, ρ inside the transition region is on the same order of magnitude as that inside the pore. Accordingly, the relevant length scale is the same as that of the pore, or . In summary, the transition region is thin due to geometrical features of the pore.

In the limit of , the transition region may present finite electric potential and charge density changes.³⁶ Therefore, in our perturbation expansion model, we approximate this region as a jump discontinuity. Defining J_left = J_Z|_Z=0⁻ and J_right = J_Z|_Z=0⁺, the charge flux in the SDL is given by where Ψ_left = Ψ_m(0⁻, τ); see eqn (2e) and (9). In contrast, the charge flux inside the pore is given by ; see eqn (2d), (18) and (19). Therefore, since the transition region is thin and thus has negligible charge storage, we ensure current conservation across it by requiring J_leftA_s = J_rightA_p, i.e.,

Eqn (23) consists of two variables, i.e., Ψ_left and Ψ_right, and we need an additional equation to solve Ψ_m.

To relate Ψ_left and Ψ_right, we require current conservation also on the left surface of the transition region (see Fig. 1b), relating its charge flux to that of the SDL. To do so, we define ρ_left = ρ_m(0⁻, τ) = 0 and ρ_right = ρ_m(0⁺, τ). The charge flux inside the transition region can be approximated by with a O(δ) error, where δ is the dimensionless thickness of the region, scaled by [small script l]_p. The current conservation relation, J_Z,int = J_left, reads

For δ ≪ 1, the right-hand side can be neglected, so ρ_right + 2Ψ_right = 2Ψ_left. Eqn (17) can then be used to deriveEqn (25) shows that when , Ψ_left = Ψ_right, which is expected for thin double layers. In contrast, when , Ψ_right = Ψ_D, which is also expected, since the potential will be close to the surface potential everywhere inside the pore. Substituting eqn (25) in eqn (23) yieldswhere is the Biot number.

3.4 Governing equation

We can now combine eqn (22) and (26) to rewrite the governing equation for centerline potential aswhere and . Here, T and ϕ are the effective dimensionless time and potential, respectively. Eqn (27a) is subjected to , and . The analytical solution of eqn (27a) yields⁴⁴where κ_n tan κ_n = Bi. Eqn (27b) can be used to derive , which can then be used to evaluate Ψ and ρ inside the pore using eqn (18) and (19). Finally, eqn (9) and (25) enable us to evaluate Ψ inside the SDL.

We emphasize that eqn (18), (19) and (27) are the key result of this paper. To the best of our knowledge, this is the first solution of the charge density and potential profiles within a cylindrical pore for arbitrary pore size in the limit of low potentials, formally justified by a regular perturbation expansion. This paper highlights that the mathematical structure of Ψ remains identical irrespective of . However, the effect of modifies the charging timescale

In the thin-double-layer limit, i.e., , , consistent with the results reported previously.² In the overlapping-double-layer limit, i.e., , , consistent with the results reported in Gupta et al.³⁶ We illustrate this timescale dependence on pore size with contour plots of the time evolution of the electric potential in Fig. 2, and rationalize its mechanism in the following discussion.

Fig. 2 Contour plots of the time evolution of the charge density in the SDL–pore region for different relative pore sizes. Charge density profiles reveal different charging timescales and screening lengths corresponding to different relative pore sizes. All the plots share the same axes.

In the thin-double-layer limit, the pore remains uncharged except in the vicinity of its surface. As shown in Fig. 2 and the ESI Video (ESI), for a_p/λ = 10, τ = 0.2 displays a charge density profile very close to that of τ = 1, indicating an earlier saturation of the profile and therefore a lower timescale. This is also backed up by the calculations, which predict τ_c ≈ 0.19 for ; see eqn (27c). In contrast, when double layers are not thin, i.e., , a radial distribution of charge propagates throughout the pore, producing a transient axial gradient of charge. The charge profiles continue to develop in both radial and axial directions and appear to saturate around τ = 1. This is also consistent with our prediction since τ_c ≈ 0.7; see eqn (27c). The charging of lies between the two other scenarios described. Finally, we also note that while the overlapping double layers take longer to charge, they also store more charge throughout the pore, resulting in a trade-off of charge density and charging timescale.

The characteristic feature that sets the charging timescale of arbitrary pore sizes is the interplay of electromigrative and diffusive fluxes. Therefore, we construct a vector plot of diffusive and electromigrative fluxes by employing our analytical derivation; see Fig. 3. First, we note that for both and , the radial diffusive flux is always balanced by the radial electromigrative flux, even though the radial diffusive and electromigrative fluxes increase with time; see Fig. 3. This is a classical feature of double-layer charging; see ref. 11 for more details. Overall, the balance in the radial direction implies that the timescale of pore charging is controlled by the fluxes in the axial direction. For , the axial gradient of charge vanishes and electromigration is the only mechanism promoting axial transport of charge, in consonance with the literature.^2,27,28,30 In contrast, for , both electromigration and diffusion cooperate in driving charge transport; see Fig. 3. However, this increase of charge flux in the overlapping double layer limit is smaller than the boost in charge density, leading to a longer charging timescale.

Fig. 3 Vector field plots of negative charge diffusive and electromigrative fluxes in the pore for different relative pore sizes. Arrow lengths are logarithmically scaled. Charge flux is only driven by electromigration for thin double layers, but set by a balance of diffusion and electromigration for overlapping double layers. All the plots share the same axes.

4 Results: potential, charge and current

In this section, we analyze the analytical predictions of the electrical potential and charge density profiles. We also validate our analytical predictions with the results from DNS.

4.1 Axial dependence of potential and charge profiles

The centerline potential, given by eqn (27b) with , and the charge density distribution, given by eqn (17), are shown in Fig. 4. We note that we obtain excellent agreement between the analytical results and DNS for both Ψ_m and ρ_m.

Fig. 4 Centerline potential and charge as a function of the axial coordinate for different . (a) and (b) for τ = 0.15. (c) and (d) for τ = 0.42. ψ_D = 0.4 and Bi = 8 for all plots. Smaller results in an increased centerline potential and charge density, and a larger potential jump at the mouth of the pore. The solid lines are predictions from the perturbation expansion model and the dashed lines are results from DNS.

Fig. 4b and d show that the centerline charge density increases monotonically as is reduced, vanishing in the thin-double-layer regime, and approaching twice the applied potential for overlapping double layers. At steady state, ϕ → 0, , and ; see eqn (17). The increase in centerline charge induces a larger potential change at the SDL–pore transition in order to balance diffusion and electromigration in this region.

4.2 Radial dependence of potential and charge profiles

The radial dependence of the pore potential and charge distribution, given by eqn (18) and (19), are shown in Fig. 5 for . The charge distribution propagates gradually along the Z-direction, being larger near the mouth of the pore. The potential and charge distributions are related in the radial direction in such a way that the net radial flux vanishes for all times. Fig. 5a and b show opposite dependences of charge density and electric potential on the axial coordinate – ρ decreases and Ψ increases with Z. Essentially, the effect of the charge flux in the SDL is to transport charge into the pore, which subsequently diffuses from the mouth to the end of pore. On the other hand, potential is screened, i.e., reduced, by the charge transported from the SDL. Thus, it is lower closer to the mouth of the pore and increases with Z. At steady state, the charge and potential distributions become independent of Z. Using eqn (18) and (19), the distributions for potential and charge throughout the pore at steady state are found to beand

Fig. 5 Potential and charge over the applied potential as a function of the radial coordinate for different axial locations. (a) and (b) For τ = 0.15. (c) and (d) For τ = 0.42. ψ_D = 0.4, a_p/λ = 2 and Bi = 8 for all plots. The solid lines are predictions from the perturbation expansion model and the dashed lines are results from DNS.

We note that smaller implies higher charge density throughout the pore. However, we reiterate that this comes at the cost of a higher charging timescale, since .

4.3 Dependence of charge flux on relative pore size

The centerline potential gives an important physical descriptor of the charging process in the dimensionless charge flux.

Eqn (2e), (25), and (27b) can be utilized to show that

The charge flux comes out to be initially independent of the relative pore size since double layers have not yet formed in the initial state, and the SDL is subjected to the potential gradient imposed by the applied potential; see eqn (3). The effect of is thus to control the charging timescale. In fact, recall that and note that J_right = J_right(T), alternatively demonstrating that narrower pores take longer to charge. This is illustrated in Fig. 6, which shows exponential-like behavior at large times, with steeper descents for larger .

Fig. 6 Dimensionless charge flux at the mouth of the pore (Z = 0⁺) vs. time with ψ_D = 0.4 and Bi = 8. Charge flux decreases monotonically as a function of relative pore size, with exponential decays at late times. The solid lines are predictions from the perturbation expansion model and the dashed lines are results from DNS.

4.4 Dependence of charge flux on the Biot number

The ratio of diffusion coefficients in the static diffusion layer and in the pore, , may not be unity. This might be due to different diffusion mechanisms in the SDL and the pore, i.e., bulk diffusion in the former vs. Knudsen diffusion in the latter due to its narrow size.⁵² The influence of this ratio is encompassed in the Biot number, , which in this case is a ratio of charge transport resistance in the pore vs. in the SDL. Fig. 7 illustrates the dependence of pore charging on the Biot number. It shows that the higher the Biot number, the lower the resistance to charge transfer in the SDL, producing an intuitive effect: enhancement of charge flux for short-times (e.g., owing to higher SDL diffusivity), but also quicker saturation of the pore charge storage. Thus, the order of the curves for different Biot numbers swaps over time. We also find that the long-time charging is given by the dominant timescale τ_c/κ₁², where κ₁ is the first eigenvalue satisfying the characteristic equation κ_n tan κ_n = Bi. This is illustrated by the good agreement between the charge fluxes and their approximation by the first mode in the Fourier series of eqn (30), represented by dash-dotted lines. It should be noted that net charge storage is not influenced by Bi, as will be shown in Section 5.2.

Fig. 7 Influence of the Biot number on the dependence of charge flux with time. Higher Biot implies faster diffusion in the SDL, with an enhancement of the initial charge flux for short-times, but also quicker saturation of the pore charge storage. Continuous lines represent analytical results and the dash-dotted lines show the first mode of the Fourier series in eqn (30). Relative pore size: .

4.5 Validity for higher potentials

For higher applied potentials, the full PNP equations assuming radial equilibrium and relaxing the constraint of constant salt follow from eqn (2) asThe main simplification that comes from the linear regime Ψ_D ≪ 1 is the neglect of the term , whence the salt transport equation resumes to a transient diffusion equation. Given the initial and boundary conditions for salt, the trivial solution s(R, Z, τ) = 2 holds, i.e., we can assume salt to be constant. For higher potentials, though, that approximation ceases to be valid and the variation of salt density in the domain influences the strength of electromigration of charge.

In order to assess the validity of our analysis, beyond which the neglect of higher-order corrections in Ψ_D becomes unwarranted, we compare our results to the DNS for different applied potentials in Fig. 8. We obtain good agreement between potentials and charge profiles until the applied potential reaches approximately twice the thermal potential, or approximately 50 mV at room temperature. Above that threshold, a non-linear potential response in the SDL is observed in the DNS, and thus the first-order corrections underpredict the potential in the pore and overpredict the charge stored. Indeed, Fig. 8e shows that the charge flux predictions agree with DNS results for early times (τ < 0.1), even for moderately high applied voltages, with an error of 8% for 200 mV (but of 15% for 100 mV). However, for late times, non-linear terms become important, introducing electromigration of salt which can affect the transport of charge.

Fig. 8 Validity of analytical results. Comparison of the proposed model (continuous lines) with DNS results (dashed lines) for different applied potentials. (a) and (b) τ = 0.07, (c) and (d) τ = 0.15. (e) Time dependence of the charge flux for a_p/λ = 2. Good agreement is obtained up to twice the thermal potential (≈50 mV at room temperature). The current shows good agreement for early times even for moderate potentials, but salt corrections to charge electromigration must be taken into account for long times.

5 Analysis of capacitance

Several studies on electrode charging invoke transmission line representations, theoretically developed for thin double layers,^27–29 to address electric potential and capacitance in these systems.^22–25 In this section, we demonstrate that a similar transmission line model can be constructed for arbitrary pore sizes. However, a time-dependent transition-region resistor makes it challenging to compare it directly with experiments. We also derive an effective capacitance for arbitrary pore sizes from the steady-state charge density profiles.

5.1 Transmission line circuit

First, we derive the value of conductivity for arbitrary . Briefly restoring dimensions and utilizing eqn (2e), (18) and (19), we find that axial charge flux J_z inside the pore is given assuch that dimensional pore conductivity is given by . Note that the conductivity for is only based on electromigrative charge flux. In contrast, the conductivity for is enhanced because both diffusive and electromigrative fluxes contribute to current; see Fig. (3).

Next, we derive the dimensional capacitance per unit surface area, C̃_p. We note that 2πa_pC̃_p(ψ_m − ψ_D) should yield the total charge stored per unit axial length. By utilizing eqn (19) and integrating along the radial direction, it is straightforward to show that . In the thin-double-layer limit , we recover the usual expression that . In contrast, in the thick-double-layer limit , we obtain . This result is intuitive since it demonstrates that the length scale of capacitance in the overlapping-double-layer limit is controlled by the pore size. Based on these expressions, we can also recover , consistent with eqn (27c).

The linear relation between average diffusive and electromigrative axial fluxes allowed us to determine expressions for pore capacitance and resistance that satisfy Ohm's law and the definition of a capacitive element. However, the key distinguishing feature from the classical thin double layer analysis^27,28 is the inclusion of potential change across the transition region as derived in eqn (25), which comes from the charge flux matching. In order to account for this change in potential, we add a resistor representing the interface, as shown in Fig. 9a. Its dimensional resistance is determined from Ohm's law,

We derive an expression for dimensionless transition-region resistance R_t(τ) in Appendix A along with the traditional expression for SDL resistance, since their final expressions are not utilized for further analysis.

Fig. 9 (a) Top panel: Transmission line circuit schematic for arbitrary pore size. (I) SDL, (II) transition region, (III) pore. SDL and transition resistances given by eqn (A2) and (A5). Bottom panel: Equivalent representation of the charge balances in an infinitesimal volume in a pore. (b) dimensionless pore capacitance per unit surface area, , (c) dimensionless charging time, . With the presence of a transition-region resistor, the transmission line model becomes more intricate and loses some applicability.

Fig. 9b and c show the pore capacitance and charging time as per the transmission line model. The sharp increase in pore capacitance per unit surface area for smaller may suggest a blowup of charge flux for overlapping double layers, which is not observed in Fig. 6. In fact, though pore capacitance increases, pore resistance decreases commensurately, in a way that yields τ_c ∼ 1 for overlapping double layers. Even more importantly, pore capacitance of arbitrary reported in our manuscript is useful for the transmission circuit analysis that predicts centerline potential, but since it's based on the potential difference Ψ_m − Ψ_D, it is not representative of experimental measures of capacitance, which are based on Ψ_D.^3,43 This is the motivation for the definition of an effective capacitance that we develop in the next section.

5.2 Macroscopic perspective

Despite the possibility of representing pore charging by the transmission line circuit described in Section 5.1, its time-dependent transition-region resistor enters as an additional factor influencing the pore charging performance. In addition, the definition of pore capacitance based on the potential difference Ψ_m − Ψ_D makes it incompatible with experimental measures. Therefore, in the interest of facilitating the analysis of pore-size effects on charging performance, our focus lies on a direct comparison via a macroscopic perspective, i.e., from characterizing the total charge stored by the pore at steady state. We employ the definition of effective volumetric capacitance C_eff as the ratio of total charge stored per applied voltage and unit volume. Going back to dimensionless variables and performing a radial integration of the steady-state charge density profile of eqn (29), it is straightforward to obtainwhere capacitance is scaled by . Eqn (35) shows that C_eff = τ_c, i.e., dimensionless volumetric capacitance is equal to dimensionless charging time; see Fig. 9c and 10. The latter shows that the volumetric capacitance decreases with an increase in relative pore size, that is, overlapping double layers present optimal energy storage. A similar behavior of energy density increase with width reduction due to a decrease of the electroneutral region has been reported in ref. 53 for nanopores, and ref. 33 for pores ranging from 0.5 to 10 nm. However, in our model, this comes at a cost of a corresponding increase in the charging time τ_c. From an engineering perspective, this prediction means that an increase in the energy density, E = C_effΨ_D², with a change in relative pore size is unaccompanied by changes in power density, P = E/τ_c; see the inset in Fig. 10. In fact, this relationship between energy and power density is consistent with experiments; see Fig. 1 of ref. 6. To the best of our knowledge, this interplay between energy and power density and its dependence on pore sizes hasn't been derived previously.

Fig. 10 Dimensionless effective volumetric pore capacitance versus relative pore size. Narrower pores account for higher volumetric capacitances. The inset shows power density versus effective volumetric pore capacitance. Gains in capacitance for a change in pore size occur at constant power density.

5.3 Effect of pore-size distribution

An important advantage of our analysis is that it is valid for an arbitrary , thus it can also be conducted for distributions of pore sizes to study the impact of polydispersity. Within the limitations of our model, i.e., specificity to non-interacting pores, we propose a simplified model to examine the influence of double-layer thickness over electrode charging. To this end, we assume a log-normal probability distribution function of pore sizes, in consonance with experimentally measured pore-size distributions,^54,55 and perform averages of the pore properties described in Section 5.2 over the pore-size distribution; see Appendix C. In our analysis, we determine the effects of pore-size average and polydispersity. Fig. 11 shows that distributions with lower averages or polydispersities of relative pore sizes present higher average capacitances (i.e., the electrode capacitance) due to elevated volumetric capacitance of pores with low relative pore size. Though this conclusion is specific to the probability density function employed, the principle of boosting average capacitance with regards to double layer thickness by increasing the frequency of narrow pores should hold in general. Optimal energy density for monodisperse pore distributions has also been reported in Monte Carlo simulations of nanopores⁵³. Nevertheless, the one caveat that also follows from our analysis is the accompanying increase in electrode charging time.

Fig. 11 (a) Dimensionless average electrode capacitance per unit volume over the electrode porosity 〈C_eff〉/ϕ vs. average relative pore size for a log-normal size distribution. (b) dimensionless averaged electrode capacitance per unit volume vs. polydispersity of the pore-size distribution, Γ, i.e., standard deviation of relative pore size over average relative pore size. An increase in the polydispersity implies a higher frequency of wider pores, which results in a decrease of the electrode capacitance.

There is an ongoing debate about the experimentally observed behavior of capacitance for sub-nanometer pores. While some studies report an anomalous areal capacitance (i.e., capacitance per unit area) increase in sub-nanometer pores,^6,56,57 other works claim that it is roughly independent of pore size in that regime,^58,59 the discrepancy being attributed to inaccuracies in BET isotherm surface area determination for subnanometer pores. The results of our work seem to support the latter hypothesis, showing only a mild increase of areal capacitance with pore size for (corresponding to a_p > 0.6 nm for 1 M electrolytes at room temperature); see Appendix B for calculations details. However, we acknowledge that our model may fail to capture intricacies of subnanometer pores, such as anomalous capacitance increases due to loss of solvation shells.^6,57 This is expected since the Poisson–Nernst–Planck equations do not take into account finite ion-radius and confinement effects,⁶⁰ which will become crucial in the subnanometer regime.

6 Conclusion

In summary, in this article:

• A regular perturbation expansion model for double layer charging at arbitrary pore sizes is proposed. The effects of arbitrary pore size include a charge flux matching condition that sets the potential change at the pore–SDL transition region;

• The proposed model predicts the potential and charge density profiles inside a nanopore. The predicted profiles using eqn (18), (19) and (27) show quantitative agreement with the results from DNS even for moderate applied potentials;

• Physical insight into the mechanisms setting capacitance and charging timescale of pores with arbitrary sizes is obtained: the influence of electromigration and charge diffusion is quantified;

• Electrical-double-layer charging for arbitrary pore sizes can be represented in the form of a transmission line circuit, but with the inclusion of a time-dependent interfacial capacitance. To mitigate this complexity, eqn (30) should be utilized to calculate charge flux;

• The electrode capacitance derived from an average of the total charge stored in the pore is able to capture some effects of pore size on pore capacitance reported in the literature.^6,59

Our methodology provides valuable insight into the effects of electromigration and diffusion in double-layer charging for arbitrary pore sizes. For thin double layers, electromigrative flux controls the charge flux and the charging timescale. In contrast, for overlapping double layers, the diffusive and electromigrative fluxes cooperate, enhancing the total charge flux. Yet, this increase in flux for overlapping double layers is less than the corresponding boost in charge density. This leads to a longer charging timescale in narrow pores, and thus a trade-off between charge stored and charging timescale. We also report a simplified model example of the influence of double layer thickness over electrode charging for non-interacting pores via phenomenological averages under a proposed log-normal distribution. In this case study, the same single-pore proportional increases of capacitance and charging timescale manifest in a distribution of isolated pores, predicting a constant electrode power density regardless of relative pore size.

Our approach can be extended to studying hybrid supercapacitors by adding reactions to the boundary conditions. The perturbation expansion analysis proposed here can also be utilized for asymmetric ionic valences¹⁹ and diffusivities,^41,42,61,62 scenarios which are commonly observed in electrochemical devices. Finally, to compare directly with cyclic voltammetry data obtained from experiments, a similar approach can also be employed for higher and/or time-dependent applied potentials.

Conflicts of interest

There are no conflicts to declare.

Appendix A: transmission line circuit resistances

In order to complete the transmission line circuit description in Fig. 9a, we derive the expressions of dimensional SDL and transition-region resistances. We start with the SDL resistance, which is given by Ohm's law asUtilizing eqn (2e) and (9), we haveFrom eqn (25), the potential difference across the transition region is given byUsing eqn (27b), it can be written asTherefore, the dimensional transition-region resistance follows from eqn (30) asWe choose the scale for the resistances, such that their dimensionless expressions readand

Fig. (12) shows a plot of the transition resistance over time for different relative pore sizes. As the relative pore size decreases, the resistance imposed by the transition region to the charge flux increases. This is due to the larger potential differences across the transition region for lower relative pore sizes, caused by their increased charge density.

Fig. 12 Resistance of the transition region as a function of time for different relative pore sizes. The resistance is larger for narrower pores and increases in time to maintain a potential difference (see eqn (25)) across the entrance region.

Appendix B: areal capacitance

Most experiments report results on a per unit surface area basis, hence we briefly present our results for areal capacitance here in order to compare them qualitatively to experiments. Denoting dimensionless areal capacitance by C_areal,eff and scaling it by , it follows from eqn (35) thatWe plot this result in Fig. 13.

Fig. 13 Dimensionless effective areal capacitance versus relative pore size. The increase in capacitance for larger pores predicted by the model developed in the current work accurately captures the trend observed in experiments for nanopores (Section III of Fig. 4 of ref. 6) and qualitatively agrees with the results of Fig. 1 of ref. 59.

Appendix C: electrode average properties

We note that the methodology presented here assumes non-interacting pores and overlooks pore-network configuration, pore intersections, connectivity, among others. To extend single-pore analysis for a distribution of non-interacting pores, we first consider a continuous log-normal distribution of pore sizes,where the parameters μ and σ are given in terms of relative pore size average and variance by and .

Next, we formally derive the averaging procedure employed in Fig. 11. The average volumetric capacitance is defined as the total charge stored inside the electrode Q divided by its total volume V and applied potential Ψ_D, compatible with what experiments denote as volumetric capacitance⁴³

Charge is the integral of the charge density over the entire electrode volume,where V_e and V_p are volumes of electrode material and pores, respectively. Furthermore, ρ = 0 inside the electrode material (an ideal conductor). Therefore, the results integral can be written as a sum over all poreswhere q_i is the charge of a pore with volume V_p,i, both known from the single-pore model. N represents the total number of pores. We insert eqn (C4) into eqn (C2) and take averages, denoted by 〈 〉, over an ensemble of sample electrode pore configurations drawn from the same pore-size distribution. In this way, we haveFor a non-interacting pore model, the average charge of the i-th pore only depends on its relative pore size, i.e.,Substituting eqn (C6) into eqn (C5), we haveNext, we multiply and divide the result by the average pore volume of the electrodewhich yieldsLastly, we write the charge of a pore in terms of its volumetric effective capacitance, q = C_effΨ_Dπa_p²[small script l]_p, change variables in the integrals from a_p to a_p/λ and define the electrode porosity ϕ = N〈V_p〉/V to getNext, we propose a definition of timescale for a distribution of pore sizes which reduces to a single-pore solution for a monodisperse distribution. To this end, note that eqn (20) can be integrated over the volume of a pore to give, using eqn (2d) and (5) for a blocking electrode, to givewhere J_right = J_Z|_Z=0⁺ and the index i denotes the i-th pore. Integrating over time and performing the change of variables T = τ/τ_c,i, we haveUsing (C11) and (C12), we find thatOur motivation in generalizing this is to construct a measure of the total charging time of the electrode, not of the individual pores. Therefore, noting that neither the initial value of the flux nor its integral on time depend on pore size – see eqn (30) – we opt to define the electrode charging timescale asaveraged over the ensemble of pore configurations, which explicitly takes the formUsing eqn (C11) and (C12) where A_p = πa_p² and changing variables from a_p to a_p/λ, we haveThrough this approach we define a physically relevant average of charging timescale that is consistent with the single-pore definition.

Acknowledgements

The authors would like to acknowledge the helpful input provided by Howard A. Stone, Daniel K. Schwartz, Adam Holewinski, Gesse Roure, Nathan Jarvey, and anonymous Referee #2. P. J. Z. would like to acknowledge the support of a project that has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 847413 and was a part of an international co-financed project founded from the program of the Minister of Science and Higher Education entitled “PMW” in the years 2020–2024; agreement No. 5005/H2020-MSCA-COFUND/2019/2. The authors acknowledge that the simulations reported in this contribution were performed using the Princeton Research Computing resources at Princeton University which is consortium of groups including the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology's Research Computing department. F. H. and A. G. would like to thank the Ryland Family Graduate Fellowship for financial support.

Notes and references

1. H. I. Becker, Low Voltage Electrolytic Capacitor, US Pat., 2800616A, United States Patent and Trademark Office, 1957.
2. P. Biesheuvel and M. Bazant, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 031502.
3. C. Zhang, K. B. Hatzell, M. Boota, B. Dyatkin, M. Beidaghi, D. Long, W. Qiao, E. C. Kumbur and Y. Gogotsi, Carbon, 2014, 77, 155–164.
4. V. Presser, C. R. Dennison, J. Campos, K. W. Knehr, E. C. Kumbur and Y. Gogotsi, Adv. Energy Mater., 2012, 2, 895–902.
5. M. Winter and R. J. Brodd, Chem. Rev., 2004, 104, 4245–4270.
6. P. Simon and Y. Gogotsi, Nat. Mater., 2008, 845–854.
7. A. Muzaffar, M. B. Ahamed, K. Deshmukh and J. Thirumalai, Renewable Sustainable Energy Rev., 2019, 101, 123–145.
8. M. Gouy, J. Phys. Theor. Appl., 1910, 9, 457–468.
9. D. L. Chapman, London, Edinburgh Dublin Philos. Mag. J. Sci., 1913, 25, 475–481.
10. O. Stern, Z. Elektrochem. Angew. Phys. Chem., 1924, 30, 508–516.
11. M. Z. Bazant, K. Thornton and A. Ajdari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 021506.
12. S. E. Feicht, A. E. Frankel and A. S. Khair, Phys. Rev. E, 2016, 94, 012601.
13. M. S. Kilic, M. Z. Bazant and A. Ajdari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 021502.
14. M. S. Kilic, M. Z. Bazant and A. Ajdari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 021503.
15. B. Balu and A. S. Khair, Soft Matter, 2018, 14, 8267–8275.
16. F. Jiménez-Ángeles, G. Odriozola and M. Lozada-Cassou, J. Chem. Phys., 2006, 124, 134902.
17. F. Jiménez-Ángeles and M. Lozada-Cassou, J. Chem. Phys., 2008, 128, 174701.
18. J. Wu, T. Jiang, D.-E. Jiang, Z. Jin and D. Henderson, Soft Matter, 2011, 7, 11222–11231.
19. A. Gupta and H. A. Stone, Langmuir, 2018, 34, 11971–11985.
20. A. Gupta, A. Govind Rajan, E. A. Carter and H. A. Stone, Phys. Rev. Lett., 2020, 125, 188004.
21. A. Gupta, A. Govind Rajan, E. A. Carter and H. A. Stone, J. Phys. Chem. C, 2020, 124, 26830–26842.
22. J. M. Black and H. A. Andreas, J. Phys. Chem. C, 2010, 114, 12030–12038.
23. M. Kaus, J. Kowal and D. U. Sauer, Electrochim. Acta, 2010, 55, 7516–7523.
24. J. Kowal, E. Avaroglu, F. Chamekh, A. Šenfelds, T. Thien, D. Wijaya and D. U. Sauer, J. Power Sources, 2011, 196, 573–579.
25. G. Madabattula and S. Kumar, J. Electrochem. Soc., 2018, 165, A636.
26. C. Lian, M. Janssen, H. Liu and R. van Roij, Phys. Rev. Lett., 2020, 124, 076001.
27. R. De Levie, Electrochim. Acta, 1963, 8, 751–780.
28. R. De Levie, Electrochim. Acta, 1964, 9, 1231–1245.
29. M. Janssen, Phys. Rev. Lett., 2021, 126, 136002.
30. P. Biesheuvel, Y. Fu and M. Z. Bazant, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 061507.
31. P. Peters, R. Van Roij, M. Z. Bazant and P. Biesheuvel, Phys. Rev. E, 2016, 93, 053108.
32. A. Levy, J. P. de Souza and M. Z. Bazant, J. Colloid Interface Sci., 2020, 579, 162–176.
33. J. Varghese, H. Wang and L. Pilon, J. Electrochem. Soc., 2011, 158, A1106.
34. S. Alizadeh and A. Mani, Langmuir, 2017, 33, 6205–6219.
35. S. Alizadeh, M. Z. Bazant and A. Mani, J. Colloid Interface Sci., 2019, 553, 451–464.
36. A. Gupta, P. J. Zuk and H. A. Stone, Phys. Rev. Lett., 2020, 125, 076001.
37. P. Wiersema, A. Loeb and J. T. G. Overbeek, J. Colloid Interface Sci., 1966, 22, 78–99.
38. R. W. O'Brien and L. R. White, J. Chem. Soc., Faraday Trans. 2, 1978, 74, 1607–1626.
39. E. Yariv, O. Schnitzer and I. Frankel, J. Fluid Mech., 2011, 685, 306–334.
40. O. Schnitzer and E. Yariv, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 021503.
41. A. S. Khair and B. Balu, J. Fluid Mech., 2020, 905, A20.
42. S. H. Amrei, G. H. Miller, K. J. Bishop and W. D. Ristenpart, Soft Matter, 2020, 16, 7052–7062.
43. M. Boota, K. Hatzell, M. Alhabeb, E. Kumbur and Y. Gogotsi, Carbon, 2015, 92, 142–149.
44. W. M. Deen, Analysis of Transport Phenomena, Oxford University Press, New York, 1998.
45. A. J. Bard and L. R. Faulkner, Electrochemical Methods: Fundamentals and Applications, Wiley, 2000.
46. J. Newman and K. E. Thomas-Alyea, Electrochemical Systems, John Wiley & Sons, 2012.
47. H. G. Weller, G. Tabor, H. Jasak and C. Fureby, Comput. Phys., 1998, 12, 620–631.
48. H. Jasak, A. Jemcov, Z. Tukovic, et al., International Workshop on Coupled Methods in Numerical Dynamics, 2007, pp. 1–20.
49. G. Chen, A. Perazzo and H. A. Stone, Phys. Rev. Lett., 2020, 124, 177801.
50. J. Lyklema, Fundamentals of Interface and Colloid Science: Soft Colloids, Elsevier, 2005, vol. 5.
51. D. F. Evans and H. Wennerström, The Colloidal Domain: Where Physics, Chemistry, Biology, and Technology Meet, Wiley-Vch, New York, 1999.
52. J. B. Rawlings and J. G. Ekerdt, Chemical Reactor Analysis and Design Fundamentals, Nob Hill Pub, LLC, 2002.
53. S. Kondrat, C. Perez, V. Presser, Y. Gogotsi and A. Kornyshev, Energy Environ. Sci., 2012, 5, 6474–6479.
54. J. Jorne and E. Roayaie, J. Electrochem. Soc., 1986, 133, 1649.
55. H.-K. Song, Y.-H. Jung, K.-H. Lee and L. H. Dao, Electrochim. Acta, 1999, 44, 3513–3519.
56. J. Chmiola, G. Yushin, Y. Gogotsi, C. Portet, P. Simon and P.-L. Taberna, Science, 2006, 313, 1760–1763.
57. J. Huang, B. G. Sumpter and V. Meunier, Angew. Chem., 2008, 120, 530–534.
58. T. A. Centeno, O. Sereda and F. Stoeckli, Phys. Chem. Chem. Phys., 2011, 13, 12403–12406.
59. A. García-Gómez, G. Moreno-Fernández, B. Lobato and T. A. Centeno, Phys. Chem. Chem. Phys., 2015, 17, 15687–15690.
60. A. A. Kornyshev, J. Phys. Chem. B, 2007, 111, 5545–5557.
61. J. Kim, S. Davidson and A. Mani, Micromachines, 2019, 10, 161.
62. B. Balu and A. S. Khair, J. Eng. Math., 2021, 129, 1–18.
63. M. Mirzadeh, F. Gibou and T. M. Squires, Phys. Rev. Lett., 2014, 113, 097701.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sm01239h

---