Charging Dynamics of Overlapping Double Layers in a Cylindrical Nanopore

The charging of electrical double layers inside a cylindrical pore has applications to supercapacitors, batteries, desalination and biosensors. The charging dynamics in the limit of thin double layers, i.e., when the double layer thickness is much smaller than the pore radius, is commonly described using an effective RC transmission line circuit. Here, we perform direct numerical simulations (DNS) of the Poisson-NernstPlanck equations to study the double layer charging for the scenario of overlapping double layers, i.e., when the double layer thickness is comparable to the pore radius. We develop an analytical model that accurately predicts the results of DNS. Also, we construct a modified effective circuit for the overlapping double layer limit, and find that the modified circuit is identical to the RC transmission line but with different values and physical interpretation of the capacitive and resistive elements. In particular, the effective surface potential is reduced, the capacitor represents a volumetric current source, and the charging timescale is weakly dependent on the ratio of the pore radius and the double layer thickness.

Porous electrodes with pore diameters between 0.5-10 nm range are commonly utilized in supercapacitors [1,2], batteries [3,4], desalination [5] and biosensors [6]. These pores are filled with ions that are attracted to the surface of the pore and form an electrical double layer. The charging dynamics of such a system is typically modeled as an RC transmission line, first proposed by de Levie [7], where the double layers (DLs) are represented as capacitive elements and the electroneutral solution is represented by a resistive element. Over the past three decades, numerous studies have been conducted that exploit some form of the RC transmission line model to understand the selfdischarge mechanism [8][9][10], the effect of pore size and shape [11][12][13][14][15][16][17][18], and the effect of electrolyte concentration [19,20], among others [21]. More recently, several studies have argued for more detailed analysis of ion transport inside the DL [4,5,[22][23][24][25][26].
The ratio of DL thickness to the pore radius impacts the ion transport inside a cylindrical pore [4,5,7,26]. The standard RC transmission line model assumes that the DL thickness is significantly smaller than the pore radius [7]. However, since the concentrated electrolytes can behave as dilute electrolytes [27] with DL thickness as large as 5-10 nm [27][28][29], and since the typical pore radii in experiments are between 0.5-10 nm [1,2], the DL thickness is often comparable to the pore size. Some reports have acknowledged the limit of the RC transmission model [4,5,30] but a model of the charging dynamics that describes the overlapping DL limit remains unavailable. In this Letter, we overcome these limitations by performing direct numerical solutions (DNS) on the Poisson-Nernst-Planck (PNP) equations. In addition, we introduce a reduced-order model that quantitatively agrees with the results from DNS and is applicable to the scenario of overlapping DLs. Specifically, we study the overlapping DL limit for small potentials and dilute electrolyte concentrations, i.e., the physical conditions when the PNP equations are valid.
We consider a cylindrical pore of radius a and length l pore where a=l pore ≪ 1. The pore is filled with a binary electrolyte such that the cation and anion valences are equal to unity [ Fig. 1(a)]. We assume that the pore wall exhibits an ideal blocking electrode condition, i.e., the normal flux of ions vanishes at the pore wall. We further assume that the diffusivities of the ions are equal and are denoted by D. We describe the cation and anion concentrations by variables c AE ðr; z; tÞ and the electrical potential by ψðr; z; tÞ, where r is the radial direction, z is the axial direction, and t is time. The mouth of the pore is at z ¼ 0. In addition, we assume that there exists a stagnant diffusion layer (SDL) of thickness l SDL and cross sectional area A SDL . The bulk is located at z ¼ −l SDL where c AE ¼ c 0 and ψ ¼ 0 [ Fig. 1(a)]. In the uncharged state, i.e., t < 0, the system is electroneutral everywhere. At t ¼ 0, the walls of the pore are set to a constant potential ψ D such that the counterions are attracted towards the pore wall to form a DL in the radial direction. During the charging process, i.e., t > 0, the charge inside the DL begins to accumulate such that the charge decreases along the length of the pore. In the fully charged state, the axial gradients vanish and the charge stored inside the DL is constant throughout the length of the pore.
We define dimensionless charge density ρ¼ðc þ −c − Þ=c 0 , salt concentration s ¼ ðc þ þ c − Þ=c 0 , time τ ¼ tD=l 2 pore and potential Ψ ¼ eψ=ðk B TÞ. The coupled PNP equations can thus be written as where the dimensionless coordinates are Z ¼ z=l pore , is the Debye length, ε is the electrical permittivity, k B is the Boltzmann constant, T is the temperature and e is the charge on an electron. We perform DNS using OpenFOAM [31]. We set Ψ D ¼ 0.4, λ=l pore ¼ 10 −3 , A SDL =ðπa 2 Þ ¼ 4, l SDL =l pore ≈ 0.5, 0.5 ≤ a=λ ≤ 20, and solve the system of Eqs. (1) [31]. The boundary condition for Ψ are (i) Ψ ¼ Ψ D at the pore walls, (ii) Ψ ¼ 0 in the reservoir, and (iii) zero normal gradient of potential at the boundaries of the SDL. The boundary condition for ρ and s are (i) zero normal flux at the pore walls, (ii) ρ ¼ 0 and s ¼ 2 at the reservoir, and (iii) zero normal flux at the boundaries of the SDL. The initial conditions are ρ ¼ 0 and s ¼ 2 everywhere with Ψ satisfying Laplace's equation with the aforementioned boundary conditions (see [31] for details).
The results from DNS indicate that the variation in potential at the center of the geometry Ψð0; Z; τÞ with Z depends significantly on a=λ; see Fig. 1(b). For the thin DL scenario, i.e., a=λ ≫ 1, the potential varies gradually across the Z ¼ 0 interface. However, for the overlapping DLs, i.e., a=λ ≤ 1, the potential changes rapidly across Z ¼ 0. In contrast, the potential near the surface of the pore ΨðR → 1; Z; τÞ increases rapidly across Z ¼ 0 for all a=λ; see Fig. 1(c). These trends suggest that the potential rapidly increases due to the presence of the electrical DL. For the thin DL, the DL is present only close to the surface and thus the rapid increase in potential across Z ¼ 0 is observed only near R → 1. On the other hand, for overlapping DLs, the DL is present throughout the pore cross section and the potential increases rapidly across Z ¼ 0 for all R.
Next, we focus on a reduced-order model inside the pore, i.e., 0 < Z ≤ 1, for the overlapping DLs scenario. In this limit, we cannot invoke electroneutrality inside the pore. Furthermore, we can not exploit the separation of length scales. However, in the linear limit Ψ D ≪ 1, we can assume that the salt concentration is constant everywhere and for all times, i.e., s ¼ 2 [31]. This observation eliminates the need to solve Eq. (1b). Moreover, since the radial flux of ρ at R ¼ 0 and R ¼ 1 vanishes, we assume the radial flux of ρ vanishes for all R. Therefore, we obtain where ρ m ðZ; τÞ and Ψ m ðZ; τÞ are, respectively, the centerline charge and potential, which are to be determined. We substitute Eq. (3) into Eq. (1c). Next, we assume that ∂ 2 =∂Z 2 ≪ ðl pore =aÞ 2 ∂ 2 =∂R 2 , and integrate with boundary where I n is the modified Bessel function of the first kind of order n. We note that as Z → 1, Eq. (4) does not capture the effects of a change in geometry. However, since a=l pore ≪ 1, we can neglect the end effects. Using Ψ ¼ Ψ m at R ¼ 0 in Eq. (4) yields a relation between centerline charge and potential: Substituting Eq. (5) in (3) and (4), it is straightforward to Finally, we substitute s ¼ 2 and Eq. (5) in (1), and evaluate at the center of the pore, to arrive at which is the governing equation for the centerline potential distribution inside the pore for overlapping DLs. Equation (7) suggests a charging timescale t c;overlap ¼ l 2 pore =(I 0 ða=λÞD), which is quite different from the corresponding thin DL limit, as we show later. The initial condition is Ψ m ðZ; 0Þ ¼ Ψ D . The boundary condition at the end of the pore Z ¼ 1 is simply ∂Ψ m =∂Z ¼ 0. Next, we focus on the boundary condition at Z ¼ 0 where the pore contacts the SDL.
We assume that the SDL is electroneutral (ρ ¼ 0), which yields a linear variation in Ψ m for Z ≤ 0; see Eq. (1c). In addition, we treat the sudden increase in the potential across Z ¼ 0 [ Fig. 1(b)] as a discontinuity. To this end, we define Ψ 1 ¼ Ψ m ð0; 0 − ; τÞ, i.e., the potential just outside the mouth of the pore, and Ψ 2 ¼ Ψ m ð0; 0 þ ; τÞ, i.e., the potential just inside the mouth of the pore. Similarly, we define ρ 1 ¼ ρ m ð0; 0 − ; τÞ ¼ 0 and ρ 2 ¼ ρ m ð0; 0 þ ; τÞ. We equate the current from the SDL with the current inside the pore at where ∂Ψ 1 =∂Z ¼ ðΨ 1 =l SDL Þl pore (recall that Ψ m is linear in SDL due to electroneutrality). To relate Ψ 1 and Ψ 2 , we invoke the condition that the charge flux in the transition region is equal to the charge flux in the diffusion layer, or where δZ is the thickness of the thin transition region.
Since δZ=l SDL ≪ 1, we get the condition ρ 2 − ρ 1 ¼ 2ðΨ 1 − Ψ 2 Þ. Physically, this implies that the charge fluxes arising from the diffusion and electromigration in the transition region balance each other to maintain current equality. By substituting ρ 1 ¼ 0 and ρ 2 ¼ ðΨ 2 − Ψ D Þ=(I 0 ða=λÞ − 1) [from Eq. (5)] in the aforementioned condition, we get By utilizing Eq. (10) in (8), we obtain where By a further change of variables ϕ ¼ Ψ m − Ψ D =I 0 ða=λÞ and T ¼ τI 0 ða=λÞ, we can rewrite Eqs. (7) and (11) as which is to be solved with boundary conditions ∂ϕ=∂Z ¼ Biϕ at Z ¼ 0 and ∂ϕ=∂Z ¼ 0 at Z ¼ 1. The analytical solution to this equation is where λ n tan λ n ¼ Bi and ϕ D ¼ Ψ D (I 0 ða=λÞ − 1)=I 0 ða=λÞ. Equation (14) provides a solution for the centerline potential Ψ m ðZ; τÞ and Eqs. (4)-(6) can subsequently be used to predict ρ m ðZ; τÞ, ρðR; Z; τÞ, and ΨðR; Z; τÞ for the overlapping DL limit. We briefly summarize the "RC transmission line" model for the thin DL limit [5,7]. This model assumes that the region outside the DL is electroneutral (ρ ¼ 0) and can be represented by an effective resistance per unit length R pore ; see Fig. 2. Consequently, the potential outside the DL is only dependent on the axial direction. The DLs are represented through an effective capacitance per unit (axial) length C pore such that the current from the capacitors in the radial direction is transported through the resistors in the axial direction; see Fig. 2. Finally, the SDL is assumed to be electroneutral and is thus also represented by an effective resistance per unit length R SDL . For the linear limit Ψ D ≪ 1, C pore ¼2πaε=λ, R pore ¼λ 2 =ðπa 2 DεÞ and R SDL ¼ λ 2 =ðA SDL D εÞ. Inside the pore, the current at an arbitrary location z can be evaluated as iðzÞ¼R pore −1 ∂ψ=∂z.

FIG. 2.
Effective circuit model is identical for the thin and overlapping DL limits. However, the values of C pore , R pore , R SDL , and ψ D;eff are different for the two limits. A summary of these values is provided in Table I. PHYSICAL REVIEW LETTERS 125, 076001 (2020)

076001-3
Furthermore, the change in axial current is equal to the current arising from the capacitor, or ∂i=∂z ¼ C pore ∂ψ=∂t.
We note that the mathematical structure of the governing equations in the two limits of DL thickness, i.e., Eqs. (13) and (15), is identical. Therefore, the effective circuit model for the overlapping DLs remains the same; see Fig. 2. However, the values of capacitance, resistances and the effective wall potential are C pore ¼ ðε=λ 2 Þπa 2 , R pore ¼ λ 2 =(πa 2 D ε I 0 ða=λÞ), R SDL ¼ 2λ 3 I 1 ða=λÞ=(aA SDL D εI 0 ða=λÞ), and ψ D;eff ¼ ψ D (I 0 ða=λÞ − 1)=I 0 ða=λÞ; see Fig. 2 and Table I. Physically, three key features differentiate the scenario of the overlapping DLs (a=λ ≲ 1) from the thin DLs (a=λ ≫ 1). First, the effective driving force for DL charging is reduced from ψ D to ψ D;eff ¼ ψ D (I 0 ða=λÞ − 1)=I 0 ða=λÞ. Second, since the DLs fill the entire pore, the capacitor represents a volumetric current source (see C pore in Table I), which is in contrast to the thin DL limit where the capacitor is a surface areal current source. Lastly, the effective resistance from the static diffusion layer is different as compared to the thin DL limit as this resistance also depends on the geometric parameter a=λ; see Table I. Overall, due to the reduced ψ D;eff , the net current predicted from the overlapping DL model would be lower than the current predicted by extrapolating the thin DL limit. Furthermore, due to the volumetric nature of capacitance, the effective time scale of charging becomes the diffusion time scale along the pore, or t c;overlap ≈ l 2 pore =D. Next, we compare the predictions of the reduced-order models with the results from the DNS simulation for a=λ ¼ 0.5 and a=λ ¼ 1 [Fig. 3]. The overlapping DL model displays excellent agreement between Ψ m ðZ; τÞ and the DNS simulations, including the rapid increase in Ψ m across Z ¼ 0; see Figs. 3(a) and 3(c). We note that the thin DL model is unable to capture the variation in Ψ m and thus should not be extrapolated in systems where a=λ ≲ Oð1Þ. Furthermore, the predicted ρ m ðZ; τÞ also exhibits excellent quantitative agreement with the DNS results; see Figs. 3(b) and 3(d). Lastly, our model is able to capture the radial variation in ρðR; Z; τÞ; see Fig. 3(e). Therefore, our model is able to capture most of the details of the DNS solution.
We plot the dimensionless current at the mouth of the pore IðτÞ, i.e., current normalized by πa 2 ec 0 D=l pore , for different a=λ [Fig. 4]. We find that IðτÞ decreases with an increase in a=λ, which implies that the utilization of pore volume for DL charging decreases with an increase in a=λ. We observe that the thin DL model is able to accurately predict IðτÞ for a=λ ¼ 5 and 20 [ Fig. 4]. However, for a=λ ¼ 0.5, the predicted values of IðτÞ from the thin DL model are significantly higher than those predicted by the DNS simulations. In contrast, our overlapping DL model is  able to quantitatively capture the trend. The reduction in the values of IðτÞ between the two models arises mainly because of the lower Ψ D;eff , an important feature not captured by the thin DL model. The analytical form of our analysis can be readily used to characterize cyclic voltammetry data in experimental investigations of batteries and supercapacitors. Furthermore, the dependence of current and charging timescales on a=λ is useful for the design of nanoporous electrodes. Our results are also informative for ion transport inside nanopores for biosensing and nanofluidic applications [38,39], and are especially relevant for ion transport in conical nanopores [40]. Our approach can be extended to account for finite ion size [41], dielectric decrement [42], electrolyte valence [43,44], electrolyte mixtures [45], asymmetric diffusivities [46] and atomistic physics [47,48]. Finally, our approach can also be extended to ion-selective electrodes and large potentials where additional effects such as concentration polarization can become important [49].