Equivalent circuit and continuum modeling of the impedance of electrolyte-filled pores

Batteries, supercapacitors, and several other electrochemical devices charge by accumulating ions in the pores of electrolyte-immersed porous electrodes. The charging of such devices has long been interpreted using equivalent circuits and the partial differential equations these give rise to. Here, we discuss the validity of the transmission line (TL) circuit and equation for modeling a single electrolyte-filled pore in contact with a reservoir of resistance $R_{r}$. The textbook derivation of the pore-reservoir impedance $R_r+Z_p$ from the TL equation does not correctly account for ionic current conservation at the pore-reservoir interface. However, correcting this shortcoming leads to the same impedance. We also show that the pore impedance $Z_p$ can be derived directly from the TL circuit, bypassing the TL equation completely. The TL circuit assumes equipotential lines in an electrolyte-filled pore to be straight, which is not the case near the pore entrance and end. To determine the importance of these regions, we numerically simulated the charging of pores of different lengths $\ell_p$ and radii $\varrho_p$ through the Poisson-Nernst-Planck equations. We find that pores with aspect ratios beyond $\ell_p/\varrho_p\gtrapprox5$ have impedances in good agreement with $Z_p$.


I. INTRODUCTION A. The physics of charging porous electrodes
Electrolyte-immersed porous electrodes are used in several technologies, including in batteries [1], solid oxide fuel cells [2], electrochemical sensors [3], supercapacitors [4,5], and capacitive deionization devices [6].In these applications, the porous electrodes typically contain pores of different shapes, widths, and lengths, connected hierarchically.When a potential difference is applied between two porous electrodes, migration of ions in electric fields leads, in each electrode, to the accumulation of one type of ion and an opposing electric charge on the electrode surface, which together are called the electric double layer (EDL) (see Fig. 1 for a schematic summary of the Introduction).Ions also diffuse if they pile up or dwindle locally and convect if the applied potential drives electroosmosis [7].Lastly, narrow pores can contain only so many finite-size ions, so the ionic fluxes are also affected by steric repulsions [8][9][10].A theoretical model for all these effects should involve at least the Poisson equation for the electrostatics, the Navier-Stokes equation for the fluid flow, and modified Nernst-Planck equations to describe the flux of finite-size ions; solvent-free ionic liquids would require a yet-to-be-developed continuum model instead [11].These equations should be solved in a porous electrode's 3d geometry, resolving charge storage in nanometre-wide pores and ionic fluxes through mesopores and between the electrodes over micrometers.Currently, computational resources do not allow one to do so.
Many models for porous electrode charging thus ignore their large-scale structure and instead focus on the charging of idealized pores, usually either a few nanometres or micrometers wide (see the second box in Fig. 1).Fluid flow is also often neglected, which is apposite for small applied potentials [12].The resulting Poisson-Nernst-Planck (PNP) equations were solved numerically [13][14][15][16][17][18] and analytically [16,17,19,20].While single-pore models oversimplify the charging of a porous electrode, numerically solving the PNP equations in a single pore is still computationally expensive, so the first mentioned studies go back less than two decades.

B. Single-pore equivalent circuit models
Historically, porous electrode charging was first studied through circuit models [21][22][23].Again, rather than an entire porous electrode, these works considered the charging of a single pore.The electrolyte in a pore has a resistance (R p ) and the electrolyte-electrode interface has a capacitance (C), but a pore does not charge like an RC circuit because the resistance and capacitance are distributed over the pore, which can be represented by cutting up R p and C connecting the pieces in the transmission line (TL) circuit-the ladder network shown in the third box in Fig. 1.In the limit of infinitely many, infinitesimally small circuit elements, the TL circuit gives rise to the TL equation [viz.Eq. ( 43)], a diffusion-type equation for the potential drop across the capacitors of the circuit.De Levie solved the TL equation for a case of a finite-length pore of constant cross-section and capacitance subject to a sinusoidal applied voltage of angular Step response z ψc(z, t) pore reservoir t Numerical: Mirzadeh [15] Yang [18] Sec.VIB Analytical: Henrique [16] Alizadeh [19] Aslyamov [20] Numerical: Keiser [37] Eloot [38] Analytical: Sec.II Sec.VIIB de Levie [23] Janssen [54] Sec.III de Levie [24] Sec.IV Posey [33] Sec.V Sec.VIC FIG. 1. Overview of modeling approaches to understand the response of porous electrodes to an applied potential.The figure mentions a few representative references by the first authors' names; see the introduction for more references and the bibliography for full information.The new contributions of this article are indicated in green.
frequency ω, yielding the pore impedance [24], where i = √ −1.The mathematical form coth (.)/ (.) is typical for diffusion in bounded geometries-it also arises for finite-length mass transfer of electroactive species to a planar electrode, where it is called the Warburg open impedance [25,26].
Equation (1) has been widely used to interpret electrochemical impedance spectroscopy (EIS) experiments on porous electrodes, often in combination with other cir-cuit elements [27][28][29][30].For instance, the impedance of a pore in contact with an electrolyte reservoir of resistance R r reads When viewing Z p and R r as circuit elements, Eq. ( 2) follows from Eq. ( 1) as the impedances of circuit elements in series can be simply added.In terms of the underlying physics, however, adding these separate pore and reservoir impedances makes less sense.De Levie's derivation of Eq. ( 1) employed a boundary condition corresponding to a counter electrode placed at the pore entrance, effectively setting the reservoir's resistance to zero.So Eq. ( 2) reintroduces the reservoir resistance after first setting it to zero.This procedure does not correctly account for ionic flux conservation at the pore-reservoir interface [see Section IV A].Still, de Levie's derivation of Eq. ( 1) is repeated unaltered in recent textbooks and reviews [26,31,32].Shortly after de Levie [23,24], Posey and Morozumi used the correct boundary condition in their study of the TL model's step response [33].One of the contributions of this article is that we show that Eq. ( 2) also follows from the TL equation using Posey and Morozumi's correct boundary condition.Figure 2 is a "complex plane plot" of Eq. ( 2), showing its real versus its imaginary part for different ω.The plot shows a 45-degree line at high frequencies, characteristic of diffusion in semi-infinite geometries, and a 90degree line at low frequencies, characteristic of capacitor charging.The transition between these two regimes occurs around the frequency ω ⋆ = π 2 /(2R p C) [34], and the other indicated formulas follow from the limits Z p (ω → ∞) = R p /(iωC) and Z p (ω → 0) = R p /3 + 1/iωC.
Equations ( 1) and ( 2) apply to a case where all the elements in the TL ladder circuit have the same resistance and capacitance; that is, the resistance and capacitance are constant along the pore.Hence, for a pore of length ℓ p , surface area A s p , and arbitrary but fixed crosssectional area A c p , in contact with a reservoir of length ℓ r and fixed cross-sectional area A c r , we have where κ is the electrolyte conductivity and c EDL is the EDL capacitance per unit electrode area.To connect Eqs. ( 1)-( 3) to the charging of an electrolyte-filled pore, R p , R r , and C must be expressed in terms of electrolyte properties and the pore and reservoir geometry.We follow the choice of most authors and consider a cylindrical pore of radius ϱ p [4,27,28,[35][36][37][38][39] and a cylindrical reservoir of radius ϱ r , so that A c p = πϱ 2 p , A c r = πϱ 2 r , and A s p = 2πϱ p ℓ p .However, we stress that the TL circuit may just as well be applied to pores and reservoirs with noncircular cross-sections.In this article, we will use the Poisson-Nernst-Planck equations to model the response of dilute electrolytes to small applied potentials.At steady state, this model yields the capacitance  c EDL = ε/λ D , where ε is the electrolyte permittivity and where the Debye length λ D is the characteristic width of the equilibrium EDL.Moreover, the PNP equations apply to electrolytes with a conductivity κ = εD/λ 2 D , with D being the ionic diffusivity, assumed to be equal among cations and anions.Inserting all these expressions into Eqs.( 1) and ( 2) seemingly gives us a theoretical impedance for arbitrary ℓ p , ℓ r , ϱ p , ϱ r , D, and λ D .This is not the case.As we explain below, underlying the derivation Eqs. ( 1) and ( 2) are several assumptions on the relation between these parameters, for instance, that the pore has a large aspect ratio (ℓ p ≫ ϱ p ) and thin EDLs (ϱ p ≫ λ D ).

C. Circuit models for porous electrode charging
Several papers extended the TL model to account for, for instance, Faradaic processes [24], contact resistances, electrodes with resistance [40], and various pore shapes [37,41].Z p and the impedances of other TL-like circuits were also connected in "super" circuits to describe the charging of porous electrodes containing different-sized [39] or hierarchically connected pores [42][43][44].Others represented porous electrodes by a parallel connection of m identical pores, for which the total impedance reads Z = R r +Z p /m [25,26,35,36].Identifying R p /m = R tot and Cm = C tot , however, yields that is, of the same form as Eq. ( 2), but with a different interpretation of its variables.
Equations ( 2) and ( 4) having the same functional form signals a general problem of interpreting EIS data by equivalent circuits: fit parameters do not always have clear interpretations.EIS on porous electrodes often yields data with shapes similar to the one in Fig. 2 [27][28][29][30][45][46][47].One can fit Eq. ( 4) to such data, for instance, with impedance.py[48] or commercial software, or one can quickly estimate and C tot ≈ π 2 /(2R tot ω ⋆ ) from the limits and the 45-to-90-degrees transition of the complex plane plot.Either way, while the complex plane plot of an electrode with thousands of pores may look like that of Eq. ( 4), unless one has verified that all assumptions underlying its derivation are satisfied, it is unclear how the fit parameters R tot , R r , and C tot relate to the microscopic details of the system at hand.By some independent experiment(s), one should thus determine the number of pores and their size and shape, verify that all pores are the same, verify that there are no hierarchical connections, etc.Until that time, the fit parameter R tot , for instance, is little more than a shorthand for 3{Re Another related problem of interpreting EIS spectra by equivalent circuits is that two circuits accounting for different mechanisms may have the same impedance.Concretely, say one studies the effect of pore shape on porous electrode charging and that a particular complex plane plot can be fitted well by the equivalent circuit model of Keiser, Beccu, and Gutjahr [37] for the impedance of different shaped pores.Such a good fit, however, does not preclude some other straight-pore model, accounting for additional physical mechanisms, from fitting the same data.

D. Microscopic models for single-pore charging
While there is a historical tradition of interpreting EIS data through equivalent circuits, the above two examples showed some of their limitations.Today's computational methods and resources allow one to predict EIS data through continuum models and molecular simulation [10,[50][51][52], which can capture hitherto neglected phenomena like image charge interaction, finite ion size, and nontrivial electrode geometries.The behavior of such complex systems might sometimes still be caught by equivalent circuits.Still, it is better to start from a first-principles model and derive its reduced-order behavior than to pose an equivalent circuit model and view its fitting to data as a justification of the model itself.
Before one can understand the EIS response of electrodes containing thousands of intricately-connected different pores, one should understand the EIS response of model geometries.In this regard, the mentioned singlepore PNP modeling studies [10,13,[15][16][17][18][19][20] helped to verify and extend the classical circuit models of de Levie and his contemporaries.In one of these works, we analytically solved the PNP equations for the charging of a single slit pore in contact with an electrolyte reservoir of negligible resistance [20].The case of small applied potentials and thin EDLs yielded an expression of the same form as the TL model's potential relaxation [viz.Eq. ( 58)] [53].In place of the TL circuit's R p C appeared ℓ 2 p λ D /(h p D), with h p the pore's width.The same expression results from multiplying the pore's capacitance C = εh p ℓ p /λ D and electrolyte resistance R = λ 2 D ℓ p /(εDh p ), both per unit length in the in-plane direction.Hence, in this case, there is an exact analytical correspondence between the microscopic 3d continuum model (PNP) and the reduced-order TL model, with an exact expression of the circuit parameters R p C in terms of electrode and electrolyte properties.That means that, in this case, the fit parameters of the TL model relate unambiguously to microscopic electrode and electrolyte properties.Other PNP modeling studies focused on the step response of pores in contact with an electrolyte reservoir [16][17][18].In these studies, the TL model predictions and the continuum data agreed decently but not precisely.
Despite these recent efforts, sixty years after de Levie's seminal papers, the charging of a single pore has still not been fully characterized.Consider again the cylindrical electrolyte-filled pore of length ℓ p and radius ϱ p filled with an electrolyte with a Debye length λ D and equal ionic diffusivities D, subject to a small sinusoidal voltage of angular frequency ω (ignore the electrolyte reservoir for now).Of this model's four length scales, ℓ p , ϱ p , λ D , and D/ω, 12 dimensionless ratios can be constructed (more will enter when an electrolyte reservoir, finite ion size, etc. are introduced).However, only three dimensionless ratios are independent; for instance, the Pecletlike parameter D/ω/ℓ p , the EDL overlap λ D /ϱ p , and the pore aspect ratio ℓ p /ϱ p .De Levie [23] implicitly discussed the product of the first two of these three ratios.For small ω, ions can keep up with the applied voltage, and EDLs are in quasi-equilibrium.For large ω, only the region near the pore mouth is charged and discharged.Accordingly, when we solve for the time-dependent potential in the pore [viz.Eq. ( 50)], we find that it varies over a frequency-dependent length ℓ ω = ℓ p / iωR p C called the penetration depth [23]; hence, the dimensionless ratio ℓ ω /ℓ p determines the extent to which the pore is charged.With Eq. ( 3) and the expressions in the lines below it, we find R p C = 2ℓ 2 p λ D /(ϱ p D) and Hence, ℓ ω /ℓ p is a product of two of the three mentioned dimensionless ratios.The EDL overlap parameter was thus already implicit in de Levie's work.Still, his results can only hold for ϱ p ≫ λ D , as overlapping EDL correspond to finite in-pore potential values at late times, which cannot be captured by the TL circuit.EDL overlap has only recently been thoroughly addressed by Henrique, Zuk, and Gupta through analytical and numerical PNP calculations [16,17].The third independent di-mensionless ratio, the pore aspect ratio ℓ p /ϱ p , has been virtually unexplored [54]-so far, most equivalent circuit and PNP studies of pore charging (implicitly or explicitly) took ℓ p /ϱ p ≫ 1, for the following reason.De Levie argued that, for the TL circuit to describe pore charging, equipotential lines in the electrolyte should be straight [23] and that short pores do not satisfy this condition (see page 372 of de Levie [24]).The second box in Fig. 1 shows equipotential lines based on numerical simulations described below (viz.Section VI).This figure shows that equipotential lines are not straight near a finite-length pore's entrance.This region will play a relatively larger role in the charging of short pores, so, indeed, the impedance of such pores cannot follow TL model predictions.

E. Overview
We comprehensively discuss single mesopore charging through ladder circuits and delineate by PNP modeling the validity of such circuits.Section II shows that the pore impedance Z p can be analytically derived from its corresponding TL circuit-we also discuss several popular TL-circuit extensions.Our derivations entirely bypass the TL-type modeling usually employed.Section III reviews two ways to go from the different ladder circuits to their corresponding TL equations.In particular, we generalize Janssen [55] to a case with Faradaic processes at the electrode surface.In Section IV, we derive the impedances of different pore-reservoir systems from their corresponding TL equations using Posey and Morozumi's pore-reservoir boundary condition.In Section V, we relate a pore's impedance to its response to a step potential.Section VI presents numerical results for the PNP equations in blocking mesopores.We determine the impedance of pore-reservoir systems with pores of different lengths and compare them to Eq. ( 2).While we focus on single pore charging, we also discuss porous electrodes in Section VII.We conclude in Section VIII.In Fig. 1, we indicate in green the locations of the new contribution of this work.We refer readers interested in practical applications of the TL model and its extensions to recent review papers [1,32,56] and textbooks [25,26,31,57].

A. Standard TL circuit
The TL model partitions the resistance R p and capacitance C of a pore into n pieces of resistance r k and capacitance c k , with k = 1, . . ., n.These elements are then connected as shown in Fig. 3.The top line in this circuit represents the pore's metallic surface, which is subjected to a small sinusoidal potential Ψ(t) = Ψ 0 sin(ωt).The bottom row represents the electrolyte in the pore and in a reservoir of resistance R r .
To determine the impedance of the circuit, we start at the last branch (n) and work our way to the reservoir resistor.The impedance of the last ladder rung reads Likewise, the impedance of the k-th rung reads The impedance of the complete circuit is then Z = R r + Z 1 ; note that Z 1 accounts for all ladder rungs.The first-order rational difference equation ( 7) previously appeared in Keiser, Beccu, and Gutjahr [37].That article considered noncylindrical pores, such that r k and c k varied along the circuit.We consider here the simpler case of a straight and homogeneous pore, for which c k = c and r k = r and thus R p = rn and C = cn.We rewrite Eqs. ( 6) and (7) with the scaled angular frequency ω = ωrc (throughout, bars indicate dimensionless quantities) and Z k = ra k /b k , with a k and b k to be determined, to The same expressions result if one takes the ratio of the top and bottom elements of the following vectors, with Equation ( 9) implies that By diagonalizing B as Bu ± = λ ± u ± , where and by using B n = P D n P −1 , where P = u + u − , we rewrite Eq. ( 10) to Hence, Z 1 = ra 1 /b 1 amounts to Next, using ω = ωR p C/n 2 , we rewrite the eigenvalues to which, inserted into Eq.( 13), yields Using that lim n→∞ (1 ± x/n) n = exp (±x), we find that is, Z p [Eq. ( 1)].The impedance of the circuit, including the reservoir resistance, then amounts to Eq. ( 2).

B. Contact resistance
To account for the resistance between a porous electrode and a current collector, we extend the TL circuit with a resistor of resistance r in the ladder's last rung, see Fig. 4. The impedance of the last rung now reads We rewrite Eq. ( 17) to with the same B as in Eq. (9b).Instead of Eq. ( 10), now which yields , and, in turn, which we denote Z con from hereon.The total resistance thus reads

C. Faradaic processes at electrode-electrolyte interface
The circuit in Fig. 5 models a pore with Faradaic (charge transfer) currents at its surface [24,29,30,58] and no dc gradients in potential and ion concentrations [59,60].The associated charge transfer resistance R F is partitioned into n pieces so that R F = r F /n. (The same circuit is used in the EIS analysis of solar cells and thin film diffusion [61]; in that context, R F is the recombination resistance.)In this case, With γ = r/r F , Eq. ( 23) reduces to Complex plane plot of a reservoir resistor connected to the pore impedance Zp [Eq.( 2), dotted] and extensions of the TL circuit accounting for a contact resistance, Zcon [Eq.(22), full line] and Faradaic processes, ZF [Eq.( 27), dashdotted].We set Rp/Rr = 10 and Rp/RF = 1.
By writing i ω′ = i ω + γ and dropping primes, we recover Eq. ( 8).Hence, Eq. ( 13) again holds, but the eigenvalues are now where we used that γ = r/r We thus find which is implicit in Eqs. ( 96), (103), and (104) of de Levie [24] and which we call the Faradaic pore impedance Z F from hereon.
The total resistance thus reads Figure 6 shows a complex plane plot of Eqs. ( 2), (22), and (27) for R p /R r = 10 and R p /R F = 1.

Ladders with large rungs
In his famous lectures, Feynman derived the impedance of an infinite LC ladder [62].Feynman ar-gued that, for large k, the impedance of successive rungs should be the same: Z k+1 = Z k .Barbero and Lelidis repeated this analysis for the TL circuit [63] with infinitely many R and C elements.Replacing r → R and c → C and setting [64], which, for small ωRC, displays Warburg-like scaling Z ∝ R/2 + R/(iωC).The crucial difference between the analyses of Barbero and Lelidis [63] and our derivation in Section II A is that we consider a pore whose overall resistance R p and capacitance C are fixed (and finite)-taking n → ∞, the resistors r = R p /n and capacitors c = C/n in our circuit become ever smaller.We can recover Barbero's result by replacing all r → R and c → C. In that case, Eq. ( 13) changes to with We write these eigenvalues in polar form, λ ± = |λ ± |e iφ± with φ ± being the arguments of the complex λ ± ; hence, From Eq. ( 29), one finds where the equality holds for ω = 0. Hence, |λ + | > |λ − | for ω > 0, which implies that, for ω > 0, in agreement with the result obtained by Feynman's method.

Distributed inductance
The derivation in Section II A allows us to study an LC network, not with finite L and C elements like Feynman did but with an overall L and C distributed over n elements, such that, again, C = cn and now also L = ln.Replacing the small resistors of Fig. 3 with small inductors of impedance iωl, we can again use Eq. ( 1) but replace R p → iωL, hence, Z = L/C coth √ −ω 2 LC, in agreement with Eq. ( 66) of Barbero and Lelidis [63].

Electrode resistance
Paasch, Micka, and Gersdorg [40] studied a transmission line with resistances in both channels, see Fig. 7. z z + dz Leaky" TL circuit for a pore with capacitive and Faradaic charging, with capacitance and resistances per unit length.
Such a circuit corresponds to a case where not only the electrolyte but also the electrode has a finite resistance, R s = nr s , with r s the small resistance of the elements in the circuit.To derive a recursion relation like Eq. ( 7) for this circuit probably requires repeated use of Y-∆ transformations.We have not yet been able to do so, so we leave this problem for future research.

Pores with varying section
Keiser, Beccu, and Gutjahr numerically solved the recursion relation Eq. ( 7) for pores with varying sections, for which c k and r k in Fig. 3 are not constant along the circuit [37].Analytically solving Eq. ( 7) for pores with varying sections will be difficult.Z k can still be written as Eq. ( 9), but i ω will depend on k.Hence, a product of k different matrices will appear, and we can no longer use B k = P D k P −1 .Progress may be possible for the particular case of a groove, for which de Levie found an analytical expression [41].

III. FROM CIRCUITS TO DIFFERENTIAL EQUATIONS
We review two ways of extracting a TL equation from its corresponding equivalent circuit.We focus on the leaky TL circuit (Fig. 5) for concreteness.z pore reservoir 9. Schematic showing the centerline potential ψc(z, t) and potential drop ψ d (z, t) in a reservoir-pore system.The TL model only models the pore region and accounts for the reservoir through a boundary condition.The centerline potential is drawn here with numerical data from Section VI.

A. De Levie's argument
De Levie's derivation [23] of the TL equation goes as follows.Figure 8 is a zoom-in of Fig. 5 without a specified start or end.Again, the top line in this circuit represents the electrode, which is at Ψ(t) everywhere.The bottom row represents the electrolyte phase, which has a centerline potential ψ c that varies along the pore.The voltage drop dψ c over a differential resistor is with R being the electrolyte resistance per unit length.For ψ c increasing in the z direction, the electric field and, hence, the ionic current point in the −z direction, explaining the minus sign in Eq. (32).For expressing the current dI that flows into the bottom line in Fig. 8, it is useful to introduce the potential drop ψ d (z, t) = Ψ(t) − ψ c (z, t) between the pore wall, which is at Ψ(t), and the center of the pore, which is at ψ c (z, t); see Fig. 9.The current that goes into a parallelconnected resistor and capacitor, with infinitesimal resistance R F / dz and capacitance C dz, respectively, then reads Rewriting Eq. ( 32) in terms of ψ d , taking a z derivative, and inserting Eq. ( 33), we find which is the TL equation for a case with homogeneous Faradaic surface conduction.Once we introduce the pore's length ℓ p , we can express the per-unit-length resistances and capacitance, R = R p /ℓ p , R F = R F ℓ p , and C = C/ℓ p .Still, the downside of the above argument is that, while it yields the correct TL equation, it does not inform on the boundary conditions that should be used.As a result, different authors solved the TL equation for different boundary conditions.Conversely, drawing a particular circuit including the first and last rungs of the ladder (like we did in Figs.3-5 and 7) fixes the boundary conditions-as we will show below, there is no room for variations.

B. Ref. [55] argument
One of us [55] showed how the TL equation, including its boundary conditions [viz.Eq. ( 48)] can be directly related to the TL circuit.The argument given there revolved around a finite-difference expression of the TL equation, including correct boundary conditions, which, in the limit n → ∞ is identical to a matrix differential equation that can also be derived directly from the TL circuit.Here, we repeat the argument for the slightly more involved circuit in Fig. 5 (and also shortly discuss the case of a circuit with a contact resistance, see Fig. 4).

Combining Ohm's and Kirchhoff 's laws for all rungs of a ladder circuit
For the circuit in Fig. 5, Ohm's law states that with Ψ(t) the potential of an external voltage source, Ψ k the potential drop over the k-th rung of the ladder, and I r k (t) the current through the k-th resistor.Kirchhoff's junction rule gives Now, the current into the k-th rung reads where Ψk (t) is the time derivative of the voltage drop across this rung.The above setup deviates from our previous work [55] in two places.First, the 1/r F term on the right-hand side in Eq. ( 37) was absent in Ref. [55] as we neglected surface conduction there.Second, the circuit in Ref. [55] contained R r rather than R r + r in the leftmost resistor.As a result, its Ohm's law corresponding to Eq. (35a) did not contain r.For consistency with Section II, we maintain this r.

TL equation for TL circuit with contact resistance
We can use the above arguments to find the corresponding equations for the circuit with a contact resistance, Fig. 4. In this case, we should omit the Ψ k (t)/r F term from Eq. ( 39) and set r F = r in Eq. (41).We find with M 3 ∈ R n×n .Similar to the above, we can show that Eq. (45b) corresponds to Using the same notation as before, in a finite difference scheme of Eq. ( 43), the boundary conditions Eqs. (46c) and (46d) reduce to ψ −1 = ψ 0 + ξ[Ψ(t) − ψ 0 ]/m and ψ m = 0, respectively.The latter condition modifies the finite difference for the second derivative at z m−1 as Combining the non-zero values ψ i in the vector ψ(t) = ψ 0 (t), . . ., ψ m−1 (t) ⊺ one can approximate Eq. ( 46) as with M 4 ∈ R m×m .After setting m = n, differences between Eqs. ( 45) and ( 47) are of subleading order in n.

IV. IMPEDANCE FROM TL EQUATIONS
Having derived TL equations from their corresponding circuits in Section III, we now derive the pore impedances Z p , Z con , and Z F from these TL equations.We highlight the differences between our derivations and those found in the literature.

A. Zp from TL equation for standard TL circuit
We start by considering a case without surface conduction (R p /R F = 0), for which the TL equation [Eq.( 43)] reduces to In the case of impedance spectroscopy with no bias potential, the wall potential reads Ψ(t) = Ψ 0 sin(ωt).By performing Laplace transformations [for a general function f (t), we write f (s) ≡ L f (t) ≡ ∞ 0 dtf (t) exp (−ts)] and using L ∂ t f (x, t) = s f (x, s)− f (x, 0), we find whose solution reads We can now find the current into the pore with Î(s) = −ℓ p ∂ z ψd (0, s)/R p [see discussion below Eq. ( 53)], giving (51) This yields the impedance Generally, the complex Laplace variable can be written as s = ς + iω.We set ς = 0 as we are interested in the steady state.Equation ( 52) is then identical to Eq. ( 2).This derivation of R r + Z p differs from the one found in the literature (both old [23] and recent [26,32]) in one crucial point: the boundary condition Eq. (48c).We showed in Section III how Eq. ( 48) is equivalent to a matrix differential equation based on combining Ohm's and Kirchhoff's laws for all the nodes of the TL circuit.Hence, the Robin boundary condition Eq. (48c) physically signals the conservation of ionic current.It is easier to see this if we rewrite Eq. (48c) in terms of the centerline potential, ψ c (z, t) = Ψ(t) − ψ d (z, t), to Here, the left-hand side gives the ionic current from the reservoir into the pore (z = 0 − ).As the TL model does not explicitly account for the reservoir at z < 0, Ohm's law for this region is expressed in terms of the total potential drop over the reservoir [ψ c (0, t)].The right-hand side of Eq. ( 53) represents the ionic current in the pore at z = 0 + .A partial derivative appears here as the ionic current in the pore is driven by an electric field −∂ z ψ c (z, t), which varies in the pore.Equation ( 53) is thus a statement of current conservation.Instead of Eq. (48c), de Levie applied a Dirichlet boundary condition at the pore-reservoir interface [23], ψ d (0, t) = Ψ(t).As this corresponds to the ξ → ∞ limit of Eq. (48c), we immediate find Î(s) = ( Ψ(s)/R p ) sR p C tanh sR p C by taking ξ → ∞ in Eq. (51).The impedance of the pore then amounts to Ẑ(s) = R p coth( sR p C)/ sR p C, which is identical to Z p .In turn, the reservoir can be reintroduced by connecting Z p in series with R r , yielding Eq. ( 52).The problem with this derivation is that the Dirichlet boundary condition fixes the local potential drop at z = 0, but one cannot enforce the potential there.Experimentally, one controls the potential difference between the pore wall and some far-away counter (and reference) electrode.Moreover, with the Dirichlet boundary condition, the physical interpretation of current conservation between the pore and the reservoir to which it is attached is lost.Interestingly, even though the usual derivation of Eq. ( 52) used wrong boundary conditions for the pore-reservoir connection, fixing this error led to the same impedance, Z p .

B. Zcon from TL equation for TL circuit with contact resistance
The circuit with a contact resistance [Fig.4] is governed by Eq. ( 46) [different from Eq. ( 48) in the boundary condition at z = ℓ p ], which is solved by instead of Eq. ( 50).Again, calculating the current by Î(s) = −ℓ p ∂ z ψd (0, s)/R p , we find the impedance in agreement with Eq. (22).
C. Faradaic pore impedance ZF from TL equation for the leaky TL circuit The TL equation of the leaky TL circuit was stated in Eq. (43).Tracing the steps we set in Section IV A, we see that Eq. (49a) changes to By writing s ′ = s + 1/(R F C), we find that Eq. ( 52) again holds, but with s replaced by s ′ , which is identical to R r + Z F [Eq. ( 27)].

V. IMPEDANCE FROM STEP RESPONSE
A system's impedance, Ẑ(s) = Ψ(s)/ Î(s), is usually measured by subjecting it to a small-amplitude sinusoidal voltage.One can also find the same impedance using any other voltage perturbation as long as 1) it contains all frequencies, and 2) the perturbation is small [66] (see also Yoo and Park [67] and Sec.3.7 of Lasia [26]).Hence, the impedance also follows from the current in response to a potential step Ψ step (t) = Ψ 0 Θ(t), with Θ(t) being the Heaviside step function, as In Section VI, we will use Eq. ( 57) to numerically determine the impedance of a continuum pore model from its step response.But first, we show how the TL model's impedance follows from its step response.
A. TL equation step response Posey and Morozumi [33] solved Eq. ( 48) for the case Ψ(t) = Ψ step (t) and found where α j with j = 1, 2, . . .are the solutions of the transcendental equation The current I(t) = ℓ p ∂ z ψ c (0, s)/R p into the pore amounts to (60) Inserting Eq. ( 60) into Eq.( 57), we find While it is not clear how Eq. ( 61) relates to R r + Z p [Eq. ( 2)], Fig. 10 shows that they overlap.This overlap can be understood for the case ξ = R p /R r → ∞, when α j = (j − 1/2)π solves Eq. ( 59), and Eq. ( 61) simplifies to Now, inserting the Weierstrass factorization of the hyperbolic cosine with complex argument [68] cosh into the right-hand side of For z = iωR p C, we then recover Z p on the left-hand side and Eq. ( 62) on the right-hand side.

B. Overlapping EDLs
The regular TL equation describes the charging of a pore whose EDLs are much thinner than the pore radius, λ D ≪ ϱ p .Henrique, Zuk, and Gupta studied the charging of pores with an arbitrary EDL thickness [16].Specifically, they analytically solved the PNP equations [viz.Eq. ( 69)] for a cylindrical pore subject to a small applied potential, for which they found the centerline potential where τ = 2ℓ 2 p λ D /(ϱ p D) × I 1 ϱ p /λ D /I 0 ϱ p /λ D and where α j with j = 1, 2, . . .are the solutions of the transcendental equation Equation (66) does not relax to ψ c (z, t) = 0 at late times when ϱ p /λ D ∼ 1.Therefore, when Henrique, Zuk, and Gupta interpreted their PNP model in terms of a TL-like ladder circuit, they had to include an interfacial resistance that grew monotonously over time [16].Moreover, they argued that the conductivity of the electrolyte in the pore changes from κ to κ H = κI 0 ϱ p /λ D /[I 0 ϱ p /λ D − 1].Hence, while the righthand side of Eq. ( 67) looks like the ratio of the pore to reservoir resistance, that interpretation only holds for ϱ p ≫ λ D , as the Bessel function factor in κ H then tends to unity.(In Ref. [16], the right-hand side of Eq. ( 67) also contained the ratio reservoir to pore diffusivities, which we consider here to be unity.)Retracing our steps of Section V A, we now find where we absorbed a factor Apart from that factor, we see that EDL overlap leads to a shift of frequencies (through τ ) as compared to Eq. ( 61).

VI. NUMERICAL STUDY OF PORE CHARGING A. Setup
We delineate the validity of the pore impedance Z p for pores of various aspect ratios by numerical simulations of their charging.As we are interested in a pore's impedance, its response to a small amplitude voltage, we ignore fluid flow, as the electroconvective term in the Navier Stokes equations is quadratic (and thus subleading) in the electric field [12].We first discuss numerical PNP simulations of pore charging in response to applied step potentials, much like Yang and coworkers [18].From these data, we determine the corresponding impedance Ẑ(s) using the method of Section V.
Figure 11 shows our system of interest.We consider two cylindrical pores of radius ϱ p and length ℓ p connected on either side of a cylindrical reservoir of radius ϱ r and length ℓ r ; all the cylinders' axes are aligned, so the whole system is axisymmetric.We use a cylindrical coordinate system r = (ρ, z, θ) with ρ and z being the radial and longitudinal coordinates, respectively.We set z = 0 at the entrance of the right pore to make comparisons to the TL model easier, as that model only explicitly treats the pore, with the effect of the reservoir captured in the boundary condition at z = 0.
We model the spatiotemporal evolution of the local electrostatic potential ψ(r, t) and the local cationic and anionic densities c ± (r, t) in our setup through the PNP FIG. 11.Schematic (not to scale) of an axisymmetric supercapacitor model consisting of two pores (left and right) connected to a reservoir (middle).Our two-dimensional numerical domain is colored grey, with the coordinate system's origin set to the right pore's entrance. equations, where ε is the electrolytic permittivity, e is the unit charge, D is the ionic diffusion coefficient (taken equal among cations and anions), and j ± are the cationic and anionic fluxes.Moreover, β = 1/(k B T ) is the inverse thermal energy, where k B is Boltzmann's constant and T is the temperature.
We consider all pore and reservoir walls blocking and set the initial ionic densities to c 0 throughout the system.At the time t = 0, we apply a potential difference 2Ψ 0 between the pores, which, due to the symmetry of our setup, is shared evenly between the pores.The following initial and boundary conditions thus apply with n being the outwards pointing normal vector at each boundary.Equation (70c) says that the respective boundaries are uncharged, which applies to dielectric materials.In Section VI C 4, we discuss a case where the boundaries Γ 1 and Γ 7 are conducting instead.
We will solve Eqs.(73a) and (73b) by the finite element method (FEM).To do so, we multiply them with test functions v and q, respectively, integrate over the domain and apply the boundary conditions Eq. ( 74).This yields their variational formulation We discretized Eqs. ( 75) and ( 76) using linear elements and solved them implicitly and coupled using a Newton solver from the FEniCS library [69].The mesh is generated using Gmsh [70], with the spatial resolution at the pore wall being 0.001ϱ p , resolving the Debye length by at least 10 grid points.The solver code and the script to generate the mesh are on this GitHub repository.We will study pore-reservoir systems with a fixed reservoir size of ℓ r /ϱ p = 20 and ϱ r /ϱ p = 10 and pore lengths of ℓ p /ϱ p = 1, 2.5, 5, 10, and 25.The ratio λ D /ϱ p represents the EDL overlap-we will consider cases for which the EDLs are thin (λ D /ϱ p = 0.01) and overlapping (λ D /ϱ p = 1).The dimensionless applied potential is set to Ψ0 = 0.1 throughout, corresponding to about 2.5 mV for systems at room temperature.73) and ( 74) for the local potential ψ(ρ, z, t) in a pore-reservoir system with λ D /ϱ p = 1/100 at t = 0.039 From the local potential ψ(ρ, z, t), we find a pore's centerline potential and potential drop, previously studied through the TL model, by ψ c (z, t) = ψ(ρ = 0, z, t) and    ψ d (z, t) = Ψ − ψ(ρ = 0, z, t), respectively.Figure 13 shows FEM solutions (lines) for ψ c (z, t) for various times, thin EDLs ϱ p /λ D = 100, and various pore lengths in the different panels.The colors in all panels refer to the same times in units of ϱ 2 p /D, where we picked colors from a purple to yellow scheme spanning the longest pore's (ℓ p /ϱ p = 25) relaxation.The shorter pores relax faster, so they are more purple.Figure 13 also shows Posey and Morozumi's TL model solution Eq. ( 58).As expected, discrepancies between both methods are most apparent for short pores.To draw Eq. ( 58), we needed to specify R p /R r for the different geometries.We approximated the pore's resistance R p and (half) the reservoir's resistance R r by

B. Step response
where κ is the electrolyte conductivity.The first term in R r is the resistance of a cylindrical resistor between two flat plates; this term is the exact resistance for cases where ϱ r = ϱ p .The second term in R r is Newman's resistance between a conducting disk and an infinitely large hemispherical electrode [71].The same resistance was later found by Hall, who identified it as the entrance resistance for ions entering a pore from a semi-infinite reservoir [72].By approximating R r by the two terms in Eq. (77b), we ensure we properly capture the reservoir resistance in the opposite limits of narrow and wide reservoirs.
Yang and coworkers [18] also studied the charging of a pore in response to a step potential through the PNP equations but did not incorporate the Newman-Hall term in R r .That article noted that Eq. ( 58) does not capture a pore's early-time charging, especially near the porereservoir interface.We found that adding the Newman-Hall resistance to R r yields better agreement between FEM solutions and Eq. ( 58), even at early times; see Fig. 13(d) and (e).Still, our expression for R r is an ad hoc combination of resistance expressions.The discrepancies that are still visible between both methods may be further reduced by using a better expression for R r and R p .Nevertheless, the impedance results discussed below [viz.Fig. 15] suggest that the TL model will never entirely capture the centerline potential's relaxation, even if one would have exact expressions for R r and R p .
Figure 14 is the same as Fig. 13 except for a different EDL overlap, ϱ p /λ D = 1.We now compare the FEM simulations of the PNP equations (lines) to Eq. (66) (dotted lines).To account for the Newman-Hall entrance resistance, we replaced the right-hand side of Eq. ( 67) with the ratio of Eqs.(77a) and (77b).Different from the case of thin EDLs, for overlapping EDLs, the late-time centerline potential transitions between −4 ∼ z/ϱ p ∼ 2 from a small value in the reservoir to a finite value in the pore.Equation (66) predicts that value to be 1/I 0 ϱ p /λ D , which amounts to 0.79 for ϱ p /λ D = 1 as considered here.
As for pores with thin EDLs, for long pores with thick EDLs, the FEM solutions and Eq. ( 66) agree decently.For shorter pores, we see that Eq. ( 66) overestimates the centerline potential.The transition region between −4 ∼ z/ϱ p ∼ 2-visible for all pores-is not resolved by Eq. (66).
C. Impedance

Numerical method
To calculate the impedance from the step response data, we modify Eq. ( 57) to where L num is a numerical realization of the Laplace transform defined by where t max is the last time of our numerical simulations.Integrating Eq. ( 79) by parts and using I = dQ/ dt, we find As t max → ∞, the first term on the right-hand side drops and Eq. ( 80) reduces to a known Laplace transform identity.
To determine Q(t) from our ψ(z, ρ, t) data, we note that Gauss's law gives access to the boundary condition between a charged conductor next to an insulator, eσ = −εn • E, with σ (m −2 ) being the surface charge number density, E = −∇ψ the local electric field, and n the normal vector into the conductor.We have n = ρ on Γ 6 , so eσ = ε∂ ρ ψ Γ6 or, in terms of the dimensionless potential, The total charge on one pore Q = e Γ6 dA σ is thus Our numerical solutions to the PNP equations give access to the dimensionless integral Q = 2π ℓp/ϱp 0 dz ∂ ρ ψ(ρ = 1, t).Putting Eqs. ( 78), (80), and (82) together, we find  2) (dotted and dashed) and PNP step voltage solutions for insulating (black lines) and conducting (green lines with circles) pore ends.The pore impedance [Eq.( 2)] is scaled to the pore resistance Rp; the PNP data is scaled to ℓp/(κπϱ 2 p ), i.e., the resistance of an isolated cylindrical electrolyte-filled pore.
where ω = ωϱ 2 p /D and Lnum {} = L num {} D/ϱ 2 p .Using that λ B = 1/(8πc 0 λ 2 D ) and that, in our PNP framework, the electrolyte's conductivity is κ = 2e 2 Dβc 0 , we find In ℓ p /(κπϱ 2 p ), we recognize the resistance R p of an ideal cylindrical pore filled with a dilute electrolyte [Eq.(77a)].Therefore, for thin EDLs, we can compare Ẑ(s)/R p = πϱ 3 p /(ℓ p λ 2 D ) Ẑnum directly to (Z p + R r )/R p [Eq. ( 2)].For thick EDLs, we will compare Ẑ(s)/R p = πϱ 3  p /(ℓ p λ 2 D ) Ẑnum to Eq. ( 68).Note that in numerically performing the Laplace transform in Ẑnum , we use Q(t) data for many more times than what we plotted in Fig. 13.Moreover, we note that the initial surface charge Q(0) in Eq. (83b) is nonzero.Physically, one applies a potential difference at t = 0 between pores by connecting them to a voltage source.The time it takes to apply this potential is set by the speed of electric signals in the external wiring.Meanwhile, the electric field in our geometry will relax accordingly on the dielectric relaxation time of the solvent, which is orders of magnitude faster than the ionic dynamics.We thus interpret Q(0) as the surface charge after the potential has been applied but before ions have moved.We determine Q(0) of the different pore-reservoir systems by a separate simulation of the Laplace equation-Eq.(73a) with its right-hand side set to zero, subject to Eqs. (74b) and (74c).2), should be Rr/Rp.We present data for Z(ωmax)/Rp, with ωmax = 10 4 , from Eq. ( 77) and a complex nonlinear least square fit of Eq. ( 2) to the numerical PNP data using impedance.py[48].

Impedance for thin EDLs
Figure 15 shows the numerically-determined impedances (black lines) for the same parameters as used in Fig. 13.This figure also shows Eq. ( 2) (black dotted lines), with R p /R r determined similarly to Section VI B. The TL model decently approximates the impedance of finite-length pores for aspect ratios beyond ℓ p /ϱ p > 5.For the smaller aspect ratios and at high frequencies, the numerical impedances deviate from the 45-degree phase angle associated with Z p , tending towards a pure capacitance (90 degrees).Notice that the high-frequency discrepancies nicely correspond to the early-time discrepancies of Fig. 13, as high frequencies in EIS correspond to fast processes.This means that improved models for R r and R p cannot fix all the TL model's problems, as changing R r /R p will merely shift Z p horizontally and not affect the high-frequency phase angle.Improved TL models should instead model the early-time nonlinear potential in the reservoir.
We compare the high-frequency limits of the numeri- cal data and analytical predictions in Table I.The second column shows Re(Z(ω max ))/R p as obtained by PNP, where we used ωmax = 10 4 , at which point Im(Z) is negligible.The third column lists the high-frequency limit of (R r + Z p )/R p , that is, R r /R p , which we determined for the respective parameters by Eq. (77).In line with our observations of Fig. 15, deviations between these two methods are larger for smaller aspect ratios.Next, we performed complex nonlinear least square fits of (R r + Z p )/R p [Eq. (2)] to the numerical PNP data using impedance.py[48], with R r /R p and R p C as fit parameters.Representative fits are shown for ℓ p /ϱ p = 1 and 2.5 with purple dashed lines.We also performed fits for all other aspect ratios, for which we list the fit parameter R r /R p in the last column of Table I.Even for the large aspect ratio ℓ p /ϱ p = 25, the numerical data and the R r /R p fit parameter differ substantially.Hence, even for the system for which the TL model was devised-a long pore subject to a small potential, in contact with an electrolyte reservoir filled with dilute electrolyte-there is no one-to-one relation between the TL model's fit parameters R r /R p and R p C on the one hand and the microscopic parameters characterizing the pore geometry and electrolyte properties on the other.Figure 16 shows the same PNP data for Z as in Fig. 15 but now scaled to Re(Z(ω → ∞)) instead of R p .This data representation corresponds more clearly to experiments on porous electrodes of various widths [4,30,42,47].Moreover, this data representation shows that decreasing the pore length leads to a smaller pore resistance; in the TL model, the pore's resistance is set by the difference between the high and low-frequency limits of the impedance, R p = 3{Re[Z(ω → 0)] − Re[Z(ω → ∞)]}.We conclude that decreasing pore length leads to impedance curves that progressively move towards that of a pure capacitor, as 1) the 45-degree line becomes shorter, and 2) the high-frequency regime deviates from 45 degrees (clearer visible in Fig. 15).

Impedance for thick EDLs
Figure 17 shows the numerically-determined impedances (black lines) and Eq. ( 68) (black dotted lines) for the same parameters as in Fig. 14.Again, the theoretical prediction performs decently for large aspect ratios but not for smaller ones.Overall, taking ϱ p /ℓ p = 5 as an example, the fit between numerics and theory is better in Fig. 15 than in Fig. 17.In our discussion of Fig. 14, we noted that EDL overlap leads to more involved centerline potentials than in the nonoverlapping case (Fig. 13): ψ c (z, t) transitions at the reservoir-pore interface from a small value in the reservoir to a finite value in the pore, even at late times.Henrique, Zuk, and Gupta's model captured the late-time in-pore centerline potential well.Conversely, the transition at the pore-reservoir interface was not captured, and the late-time centerline potential of short pores was overestimated.These two points may have led to the larger discrepancies between numerics and theory for short pores in Fig. 17 than in Fig. 15.

Impedance of pores with conducting ends
So far, we discussed pores whose cylindrical surface was conducting but whose ends (Γ 1 and Γ 7 in Fig. 11) were insulating.That boundary condition corresponds to the experiments of Eloot and coworkers [73] on pores drilled into stainless steel and insulating plexiglass at their ends.Conversely, pores in supercapacitor have conducting carbon surfaces on all sides except their opening.To describe such pores, we change Eqs.(70b) and (70c) to As a result, Eqs. (74b) and (74c) change to  68) is scaled to the pore resistance Rp [into which we absorbed an EDL overlap dependent prefactor, see below Eq. ( 68)], and the PNP data is scaled to ℓp/(κπϱ 2 p ), i.e., the resistance of a cylindrical pore with an electrolyte at infinite dilution.

Equation (75) changes to
Equation ( 81) changes to and Eq. ( 82) for the total charge on one pore Q = e Γ6,Γ7 dA σ becomes Figures 15 and 17 show the impedance of pores with conducting ends of various lengths (green lines with circles) as obtained from PNP solutions.For ℓ p /ϱ p = 10 and 25, these data hardly differ from the impedance of pores with insulating ends.For shorter pores, differences between both boundary conditions appear, which makes sense as a relatively larger part of the pore's charged surface area comes from its end.For short pores and thin EDLs [Fig.15], the data differ mainly at low frequencies; for thick EDLs [Fig.17], they differ mainly at high frequencies.

A. The term "Diffusion impedance"
In the context of pore charging through EDL formation at blocking electrodes, the commonly-used terminology "diffusion impedance" is a misnomer [56].As we showed in Fig. 12 (see also Fig. (3) of Henrique, Zuk, and Gupta [16] [74]), ions flow into a pore by electromigration, not diffusion.Dropping the all-important electromigration terms in the PNP equations yields a regular ionic diffusion equation.Hence, solving the ionic diffusion equation to find an electrode's impedance, as was done, for example, in Ref. [75], does not account for the relevant physics (as these authors acknowledged).Nevertheless, Ref. [75] found sensible impedances from the perturbed ion densities.How can this be?We have seen in this article that de Levie's transmission line model, a diffusion-type equation for the potential drop ψ d , accurately describes the relaxation of a pore's centerline potential.Solving a diffusion equation for ionic species and determining the impedance from the perturbed densities yields the correct impedance, as the mathematical form of all the equations is the same as the ones we used to derive the pore impedance Z p from the TL equation (but in 3d).Hence, Ref. [75] solved the correct diffusion-type equation, but the diffusing quantity is the centerline potential ψ d , not the ions.So far, we have discussed charging a single cylindrical pore in contact with a large reservoir.Different models were proposed to go from known single-pore charging behavior to predict the charging of a complete porous electrode.Here, we compare two models for an electrode with m pores.
Several papers treated porous electrodes as a bundle of m cylindrical pores connected in parallel [26,35,36]; see Fig. 18(a).In this case, the impedance of both electrodes and reservoir amounts to The current in response to a step potential, for which . The relaxation time of this system is set by the zeros of Z tot , that is, by the solution to Substituting which, up to the factor m, is the same as in Janssen [55] [and Eq. ( 59) here].An approximate solution based on Padé approximation reads α −1 j ≈ 1/3 + mR r /(2R p ), which yields the relaxation time With R p = ℓ p /(κπϱ 2 p ), C = ε2πϱ p ℓ p /λ D , and κ = εD/λ 2 D we find R p C = 2λ D ℓ 2 p /(Dϱ p ).To express the reservoir resistance R r = ℓ r /(κA c ), we equate the reservoir's cross-sectional area A c (perpendicular to the pores) to that of the pore-bundle electrode.Assuming no space to be left between the pores, each having a radius ϱ p , yields A c = mπϱ 2 p .Collecting terms, we find which, notably, does not depend on m.Lian and coworkers [76] recently proposed an alternative model for porous electrode charging.In their "stack electrode" model [Fig.18(b)], the two porous electrodes of a supercapacitor, separated by 2L and both of width H, are represented by m flat electrode "sheets" spaced h apart [so that H = h(m−1)] [76][77][78].Of these sheets, the outer ones are blocking, while the others are fully permeable to ions.Upon applying a potential difference to the two porous electrodes, with each sheet in an electrode at the same potential, ions move perpendicular to the sheets and through them, forming EDLs on both sides of each sheet (except the outer sheets).When the lateral size of the sheets is much larger than the width 2H + 2L of the setup, the potential and ion densities depend only on the coordinate z perpendicular to the sheets.Lian and coworkers [76] solved the PNP equation in this effectively one-dimensional geometry to determine each sheet's time-dependent surface charge.They showed that the stack electrode model relaxes, for small applied potential, with the same timescale as a discrete TL circuit (Fig. 3) with m rungs, with r i = r and c i = 2c for i = 1, . . ., m − 1 and c m = c, and total resistance R = i r i and capacitance C = i c i .For m ≫ 1, this circuit relaxes with almost the same timescale as the regular finite-m TL circuit, whose timescale reads [55] τ ≈ Comparing the two models in Fig. 18 and identifying L → ℓ r /2, H → ℓ p , and h → ϱ p , Eq. (94) becomes τ = 2(m − 1) obviously, with differences to Eq. ( 96) being subleading in m.For a stack electrode model whose last plate is permeable as well, both models have identical charging times.
The m parallel pores and stack electrode models both utilize the TL circuit, but they do so differently.The m parallel pores model uses Z p , which we found from the TL circuit in the n → ∞ limit.In other words, the m parallel pore model uses the n → ∞ circuit m times.By contrast, m is kept finite in the stack electrode model, with no corresponding n → ∞ limit.
While the relaxation times of both models are thus the same, their impedances are not, as we saw by comparing Eq. ( 90) to R r + 2Z 1 , with Z 1 from Eq. ( 13).This is unsurprising as the parameter m plays different roles in both models.In the m parallel pore model, increasing m corresponds to using electrodes with a larger crosssectional area.The stack electrode model, by contrast, is one dimensional, so it models a porous electrode per unit cross-sectional area.Increasing m in the stack electrode model corresponds to using thicker electrodes (if the pore width h is kept fixed) or using narrower pores (if the electrode thickness H is kept fixed).

VIII. CONCLUSIONS
We derived the pore impedance Z p directly from its corresponding TL circuit-to our knowledge, sidestepping the TL equation or other diffusion-type PDEs for the first time.As the TL circuit and its extension find use in interpreting various electrochemical devices such as batteries and fuel cells [1,2,79], our methods could be useful more broadly than for the example of EDL capacitors with porous that we focussed on here.Future work could generalize our calculations to determine the impedance of a groove [41], an arbitrarilyshaped pore [37], or to find the impedance of a case with finite electrode resistance [40].
There are at least four lengthscales relevant to the charging of a cylindrical pore: its length ℓ p and radius ϱ p , the width λ D of the EDL, and the combination D/ω of the ionic diffusion constant to the angular frequency of the harmonic voltage source.Two of the three independent dimensionless combinations of these lengthscales had been characterized.De Levie showed that a dimensionless penetration depth ∝ Dϱ p /(ωℓ 2 p λ D ) sets the characteristic length until ionic density profiles in a pore are perturbed [23]; Henrique, Zuk, and Gupta studied the effect of the EDL overlap ϱ p /λ D on pore charging.This left one dimensionless ratio, the pore aspect ratio ℓ p /ϱ p , which had received little attention.Accordingly, we studied the charging of pores of various aspect ratios by numerical simulations of the Poisson-Nernst-Planck (PNP) equations.We found impedances of long pores to agree well with Z p .By contrast, deviations were visible at high frequencies for pores with aspect ratios less than ℓ p /ϱ p = 5.Our findings are thus in qualitative agreement with Eloot and coworkers [73], who found that their experimental pore impedance data could not be fitted by equivalent circuits when ℓ p /ϱ p < 2.
The shapes of the impedance curves that we found are not unique to short pores; similar curves resulted, for instance, from an equivalent circuit model accounting for the outer surface of a porous electrode through a parallel connection of Z p and another capacitor [28].Figure 20 of that article contains experimental impedance data for a porous gold electrode; the shape of their impedance is very similar to ours in Fig. 15 for ℓ p /ϱ p = 2.5.Hence, above-45 degrees high-frequency phase angles may be explained by at least two distinct phenomena: pore aspect ratio or outer surface capacitance.Deciding which applies would require further impedance spectroscopy on different electrodes or different experiments.
We see the following directions for future work.First, an outstanding challenge is to analytically solve the PNP equations we solved numerically in Section VI.In previous work, we analytically solved the PNP equation for a long pore and negligible reservoir resistance [20].Relaxing these restrictions to describe a short pore next to a nonnegligible reservoir will be challenging.Second, the boundary conditions of the PNP equations can be adapted to pores with curved [37], rough [41,[80][81][82][83][84], or nonblocking surfaces [85,86].Third, the PNP model should be extended with finite ion sizes and dispersion and image charge interactions when pores are very narrow [10,87] or large potentials are applied, for instance, when probing a system's nonlinear impedance [88,89] or its impedance around a large bias voltage.Large applied potentials cause diffusive salt transport not captured by the TL model [20], so a pore's impedance will deviate from Z p .Last, this article aimed at bringing equivalent circuit and continuum modeling of electrolyte-filled pores closer together.It would be interesting to do the same for the equivalent circuit models and molecular dynamics simulations [90][91][92]; that is, to pinpoint the meaning of fit parameters when the TL model is fitted to molecular dynamics data.
FIG. 12. Snapshot of numerical PNP solutions (a,c) and the ionic flux contributions (b,d) for a pore-reservoir system for ℓr = 20ϱp, ϱr = 10ϱp, ℓp = 5ϱp and λD/ϱp = 1/100 at t = 0.039 (a,b) and λD/ϱp = 1 at t = 52 (c,d).In panel (b), the diffusive flux vector is scaled 100 times larger than the electromigrative flux vector.Only one pore and the immediate vicinity of the reservoir are shown here.

Figure 12 (
Figure 12(a) and (c) show numerical solutions to Eqs. (73) and(74) for the local potential ψ(ρ, z, t) in a pore-reservoir system with λ D /ϱ p = 1/100 at t = 0.039 (a) and with λ D /ϱ p = 1 at t = 52 (c).The figure shows that isopotential lines are not parallel to the pore's wall close to its entrance and end.At this intermediate time, the pore has attracted counterions and developed EDLs near the reservoir-pore interface.Figure 12(b) and (d) correspond to the same parameters as panels (a) and (c), respectively, and show the diffusive (red arrows) and electromigrative (blue arrows) contribution to the ionic fluxes, that is, the first and second terms on the righthand side of jρ,+ − jρ,− = −∂ ρ(c + − c+ ) − (c + + c− )∂ ρ ψ.Note that only in panel (b), corresponding to the same early time as in panel (a), we stretched the red arrows a hundredfold to make them visible compared to the blue arrows.Hence, the ionic fluxes are almost entirely caused by the electric field, not by diffusion.In panel (d), corresponding to the same parameters as panel (c) (overlapping EDLs and a late time), the strongest diffusive and electromigrative fluxes are near the pore-reservoir interface, where they nearly balance each other.

4 − 25 FIG. 15 .
FIG.15.Impedance Z of pores of different aspect ratios ℓp/ϱp with thin EDLs (λD/ϱp = 0.01) and other parameters as in Fig.13.The data corresponds to Eq. (2) (dotted and dashed) and PNP step voltage solutions for insulating (black lines) and conducting (green lines with circles) pore ends.The pore impedance [Eq.(2)] is scaled to the pore resistance Rp; the PNP data is scaled to ℓp/(κπϱ 2 p ), i.e., the resistance of an isolated cylindrical electrolyte-filled pore.

FIG. 16 .
FIG.16.Impedances Z of pores of different lengths, determined from PNP step voltage solutions.This is based on the same data as the black lines in Fig.15, but now all curves are scaled to their high-frequency limit, Re(Z∞) = Re(Z(ωmax)).

B
. Towards porous electrodes: m parallel pores vs. stack electrode model