Charging Dynamics of Electrical Double Layers Inside a Cylindrical Pore: Predicting the Effects of Arbitrary Pore Size

Porous electrodes are found in energy storage devices such as supercapacitors and pseudocapacitors. However, the effect of electrode-pore-size distribution over 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.


I. INTRODUCTION
Batteries and fuel cells are traditional porous-materialbased electrochemical devices. Over the past 60 years, new devices such as supercapacitors [1] have started to emerge. They are comprised of porous electrodes -typically made up of dispersions of activated carbon spheresimmersed 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 [3]. 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 [4] and memory protection in electronic devices [5]. 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 * Corresponding author: ankur.gupta@colorado.edu double layers, and another metal oxide pseudocapacitorlike which accumulates energy by performing fast oxidation surface reactions [6]. 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 [6].
The earliest models for electrode charging address the geometry of an electrolyte between flat plates, dating back to the works of Gouy [8], with later contributions from Chapman [9] and Stern [10] 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. [11] 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. [12] 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 [15], and ion correlations [16][17][18] have been studied, for instance in some continuum models [19][20][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][23][24][25][26], stemming from the pioneering work of de Levie [27,28]. Later, Biesheuvel and Bazant [2] 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] 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 [33], while others propose transient transport-equationbased 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 [36], 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][38][39][40][41][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.

II. 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 ℓ p and radius a p , where the radial and axial directions are denoted by r and z, respectively. The SDL, of length ℓ 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, −ℓ s , t) = c 0 and ψ(r, −ℓ s , t) = 0. At t = 0, the concentration of ions everywhere is c ± (r, z, 0) = c 0 . 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 ap ℓp ≪ 1, ℓs ℓp = O(1), 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., Dp Ds -dictates the interaction between the pore and the SDL. We denote Debye length by λ such that λ = εkB T 2e 2 c0 , 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 ap λ . 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][45][46] where N ± are the ion fluxes. Inside the pore and outside the pore We note that N ± = e r N ±r + e z N ±z by angular symmetry. We introduce dimensionless charge density ρ = 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 ℓs/ℓp and 1, respectively. The radii of SDL and the pore are as/ap and 1, respectively, and the Debye length is λ/ap.
c0 , salt concentration s = c++c− c0 , potential Ψ = eψ kB T , time τ = tDp ℓ 2 p , axial position Z = z ℓp , radial position R = r ap , and gradient operator∇ = ℓ p ∇. Note that − ℓs ℓp ≤ Z ≤ 1. In addition, 0 ≤ R ≤ 1 for 0 ≤ Z ≤ 1 (i.e., the pore region), and 0 ≤ R ≤ as ap for − ℓs ℓp ≤ Z < 0 (i.e., the static diffusion layer region; see Fig. 1b). With these substitutions, the set of Eqs. (1) becomes and outside the pore At τ = 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, Eqs. (2) are subjected to the following initial conditions The boundary conditions (BCs) at Z = − ℓs ℓp 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 We note that Eq. (7c) represents the ideally conducting electrode BC. Finally, at the boundaries of the SDL, i.e., R = as ap and Z < 0, the BCs are vanishing gradients, or This 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 Eqs. (2) -(8) numerically using Open-FOAM [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 [44] to obtain an analytical expression for ρ(R, Z, τ ) and Ψ(R, Z, τ ).

III. 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.

B. 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., We also introduce variables ρ m (Z, τ ) and Ψ m (Z, τ ) that represent the centerline charge and density. As per our proposed perturbation expansion, In this article, we only focus on the leading-order and first-order terms.
Inserting the expansions provided in Eq. (10) (with ρ 0 = 0, s = 2 and Ψ 0 = 0) in Eq. (2a) and (2d) and collecting the first-order terms in Ψ D yields Similarly, Poisson's equation -Eq. (1b) -after substitution of the perturbation expansions takes the form Only the known leading-order term in the salt concentration enters Eq. (11), such that Eqs. (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 Eq. (2b). Since ap ℓp ≪ 1, the diffusion in the radial direction is much faster than in the axial direction and we can assume quasiequilibrium in the radial direction [34]. Thus, with symmetry and ideally blocking conditions (see Eqs. (6) -(7)), we have J R = 0. Utilizing this condition in Eq. (2d) and collecting the first-order terms in the perturbation expansion, we get which implies that [36] where ρ m1 (Z, τ ) and Ψ m1 (Z, τ ) are to be determined. Next, by utilizing Integrating Eq. (15) with ∂Ψ1 ∂R R=0 = 0 and Ψ 1 (Z, 1, τ ) = 1 yields where I n is the modified Bessel function of the first kind of order n. By substituting R = 0 in Eq. (16), we obtain Next, by combining Eqs. (16) - (17), substituting the leading-and first-order coefficients in Eq. (10), and neglecting higher-order terms, we get Similarly, substituting Eqs. (16) and (17) into Eq. (14) and neglecting higher-order terms in the perturbation expansion of ρ, we get Next, we note that since J R = 0, so we can simplify Eq. (11) to obtain To determine the axial dependence, we average Eq. (20) over the cross-sectional area of the pore (integrating across the radial direction) by utilizing the known radial dependence in Eqs. (18) and (19). Mathematically, We emphasize that this crucial step enables us to derive a solution for arbitrary ap λ , thereby bridging the previously reported trends in thin and overlapping double layer limits [2,27,28,36]. Now, substituting Eqs. (18) and (19) into Eq. (21), we obtain The boundary conditions for Eq. (22) include Ψ right = Ψ m (0 + , τ ) and ∂Ψm ∂Z Z=1 = 0, where Ψ right , i.e., the centerline potential to the right of the SDL-pore interface, is to be determined.

C. 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 ℓ p ), we emphasize that charge density is related to the length scale of the charge gradients by the Poisson equation; see Eq. (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 ℓs ℓp = O(1). Therefore, using Eq. (2c) and as- , 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 ap ℓp , which implies ρ ΨD = O(1) [36,44,[49][50][51]; see also Eq. (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 δ = O ap ℓp ≪ 1. In summary, the transition region is thin due to geometrical features of the pore.
In the limit of ap λ = O(1), the transition region may present finite electric potential and charge density changes [36]. 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 (2e) and (9). In contrast, the charge flux inside the pore is given by (18) and (19). Therefore, since the transition region is thin and thus has negligible charge storage, we ensure current conservation across it by requiring Eq. (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 J Z,int = −( where δ is the dimensionless thickness of the region, scaled by ℓ 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 . Eq. (17) can then be used to derive Eq. (25) shows that when ap λ ≫ 1, Ψ left = Ψ right , which is expected for thin double layers. In contrast, when ap λ ≪ 1, Ψ right = Ψ D , which is also expected, since the potential will be close to the surface potential everywhere inside the pore. Substituting Eq. (25) in Eq. (23) yields where Bi = Ds Dp ℓp ℓs is the Biot number.

D. Governing Equation
We can now combine Eqs. (22) and (26) to rewrite the governing equation for centerline potential as . Here, T and φ are the effective dimensionless time and potential, respectively. Eq , which can then be used to evaluate Ψ and ρ inside the pore using Eqs. (18) and (19). Finally, Eqs. (9) and (25) enable us to evaluate Ψ inside the SDL.
We emphasize that Eqs. (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 ap λ . However, the effect of ap λ modifies the charging timescale In the thin-double-layer limit, i.e., λ ap ≪ 1, t c ≈ 2λ ap ℓ 2 p Dp , consistent with the results reported previously [2]. In the overlapping-double-layer limit, i.e., λ ap ≫ 1, t c ≈ ℓ 2 p Dp , consistent with the results reported in Gupta et. al. [36]. 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.
In the thin-double-layer limit, the pore remains uncharged except in the vicinity of its surface. As shown in Fig. 2 and the Supplementary Video, 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 ap λ = 10; see Eq. (27c). In contrast, when double layers are not thin, i.e., ap λ = 2, a radial distribution charge forms 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 Eq. (27c). The charging of ap λ = 5 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 ap λ = 2 and ap λ = 10, 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 ap λ = 10, 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 ap λ = 2, both electromigration and diffusion cooperate in driving charge transport; see Fig. 3. However, this increase in charge flux of in the overlapping double layer limit is smaller than the boost in charge density, leading to a longer charging timescale.  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.

A. Axial Dependence of Potential and Charge Profiles
The centerline potential, given by Eq. (27b) with Ψ m = φ + ΨD I0( ap λ ) , and the charge density distribution, given by Eq. (17), are shown in Fig. 4. We note that we obtain excellent agreement between the analytical results and DNS for both Ψ m and ρ m .
Figs. 4b and 4d show that the centerline charge density increases monotonically as ap λ is reduced, vanishing in the thin-double-layer regime, and approaching twice the applied potential for overlapping double layers. At steady state, φ → 0, Ψ m → ΨD I0( ap λ ) , and ρ m → − 2ΨD I0( ap λ ) ; see Eq. (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.

B. Radial Dependence of Potential and Charge Profiles
The radial dependence of the pore potential and charge distribution, given by Eqs. (18) and (19), are shown in Fig. 5 for ap λ = 2. 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. Figs. 5a and 5b 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 Eqs. (18) and (19), the distributions for potential and charge throughout the pore at steady state are found to be and We note that smaller ap λ implies higher charge density throughout the pore. However, we reiterate that this comes at the cost of a higher charging timescale, since

C. 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.
Eqs. (2e), (25), and (27b) can be utilized to show that sin 2κ n 2κ n + sin 2κ n exp(−κ 2 n T ). (30) 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 Eq. (3). The effect of ap λ is thus to control the charging timescale. In fact, recall that T = τ 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

D. Dependence of Charge Flux on the Biot Number
The ratio of diffusion coefficients in the static diffusion layer and in the pore, Ds Dp , 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 [52]. The influence of this ratio is encompassed in the Biot number, Bi = AsDsℓp ApDpℓs , 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 swap over time. We also find that the long-time charging is given by the dominant timescale τ c /κ 2 1 , where κ 1 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 Eq. (30), represented by dash-dotted lines. It should be noted that net charge storage is not influenced by Bi, as will be shown in Sec. 5.2.

E. Validity for Higher Potentials
For higher applied potentials, the full PNP equations assuming radial equilibrium and relaxing the constraint of constant salt follow from Eqs. (2) as The main simplification that comes from the linear regime Ψ D ≪ 1 is the neglect of the term ρ ∂Ψ ∂Z , 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. Good agreement is obtained up to twice the thermal potential (≈ 50mV 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.

V. ANALYSIS OF CAPACITANCE
Several studies on electrode charging invoke transmission line representations, theoretically developed for thin double layers [27][28][29], to address electric potential and capacitance in these systems [22][23][24][25]. In this section, we demonstrate that a similar transmission line model can be constructed for arbitrary pore sizes. However, a timedependent 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.

A. Transmission Line Circuit
First, we derive the value of conductivity for arbitrary ap λ . Briefly restoring dimensions and utilizing Eqs. (2e), (18) and (19), we find that axial charge flux J z inside the pore is given as such that dimensional pore conductivity is given bỹ Note that the conductivity for ap λ ≫ 1 is only based on electromigrative charge flux. In contrast, the conductivity for ap λ ≪ 1 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 pCp (ψ m − ψ D ) should yield the total charge stored per unit axial length. By utilizing Eq. (19) and integrating along the radial direc- In the thin-double-layer limit ap λ ≫ 1, we recover the usual expression thatC p → ε λ . In contrast, in the thickdouble-layer limit ap λ ≫ 1, we obtainC p → 2ε ap . 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 t c = 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 Eq. (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 transitionregion resistance R t (τ ) in Appendix A along with the traditional expression for SDL resistance, since their final expressions are not utilized for further analysis.
Figs. 9b and 9c 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 ap λ 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 ap λ 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.

B. Macroscopic Perspective
Despite the possibility of representing pore charging by the transmission line circuit described in Sec. V A, 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 Eq. (29), it is straightforward to obtain where capacitance is scaled by ε λ 2 . Eq. (35) shows that C eff = τ c , i.e., dimensionless volumetric capacitance is equal to dimensionless charging time; see Figs. 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 Refs. [53] for nanopores, and [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 Ψ 2 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.

C. Effect of Pore-Size Distribution
An important advantage of our analysis is that it is valid for an arbitrary ap λ , 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 Sec. V B over the poresize 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 [53]. Nevertheless, the one caveat that also follows from our analysis is the accompanying increase in electrode charging time.
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 ap λ > 2 (corresponding to a p > 0.6 nm for 1 M electrolytes at room temperature); see Appendix A 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 [60], which will become crucial in the subnanometer regime.  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.

VI. 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; a) b) FIG. 11. a) Dimensionless average electrode capacitance per unit volume over the electrode porosity C eff /φ vs. average relative pore size ap λ 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.
• The proposed model predicts the potential and charge density profiles inside a nanopore. The predicted profiles using Eqs. (18), (19) and (27) show quantitative agreement with the results from DNS even for moderate applied potentials; • Physical insight into the mechanisms setting capac-itance 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 timedependent interfacial capacitance. To mitigate this complexity, Eq. (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 tradeoff between charge stored and charging timescale. We also report a simplified model example of the influence of double layer thickness over electrode charging for noninteracting 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 [19] 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. 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 Sk lodowska-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.
(A5) We choose the scale ℓpλ 2 εDpAp for the resistances, such that their dimensionless expressions read and . (A7) Fig. (A1) 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. A1. 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 Eq. (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 Eq. (35) that C areal, eff = I 1 ( . (B1) We plot this result in Fig. (A2).  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 lognormal distribution of pore sizes, where the parameters µ and σ are given in terms of relative pore size average and variance by ap λ = exp(µ+ σ 2 2 ) and Var( ap λ ) = [exp(σ 2 ) − 1] exp(2µ + σ 2 ). 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 [43] 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 pores where 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 Eq. (C4) into Eq. (C2) and take averages, denoted by , over an ensemble of sample electrode pore configurations drawn from the same poresize distribution. In this way, we have For a non-interacting pore model, the average charge of the i-th pore only depends on its relative pore size, i.e., Substituting Eq. (C6) into Eq. (C5), we have Next, we multiply and divide the result by the average pore volume of the electrode which yields Lastly, we write the charge of a pore in terms of its volumetric effective capacitance, q = C eff Ψ D πa 2 p ℓ p , change variables in the integrals from a p to a p /λ and define the electrode porosity φ = N V p /V to get (C10) Next, 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 Eq. (20) can be integrated over the volume of a pore to give, using Eq. (2d) and Eqs. (5) for a blocking electrode, to give where 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 have Using (C11) and (C12), we find that τ c,i = J right,i (τ = 0) ∞ 0 J right,i (T (τ )) dT