Curvature affects electrolyte relaxation: studies of spherical and cylindrical electrodes

TL;DR

This paper investigates how electrode curvature influences electrolyte relaxation times using minimal models for spherical and cylindrical electrodes, revealing that relaxation times depend on both electrode radii rather than just their difference.

Contribution

It extends existing models to curved electrodes, demonstrating that ionic relaxation times are affected by electrode radii, not only electrode separation.

Findings

Relaxation times depend on both electrode radii.

Equivalent circuit models accurately predict ionic relaxation times.

Curvature significantly influences electrolyte response.

Abstract

With two minimal models, I study how electrode curvature affects the response of electrolytes to applied electrostatic potentials. For flat electrodes, Bazant et al. [Phys. Rev. E. 70, 021506 (2004)] popularized the "RC" timescale $\lambda_{\textrm{D}} L/D$, with $\lambda_{\textrm{D}}$ being the Debye length, $2L$ the electrode separation, and $D$ the ionic diffusivity. For thin electric double layers near concentric spherical and coaxial cylindrical electrodes, I show here that equivalent circuit models again predict the correct ionic relaxation timescales. Importantly, these timescales explicitly depend on both electrode radii, not simply on their difference.

Equations

$$C^{\mathrm{s}}=\frac{4\pi \epsilon }{{\lambda }_{\mathrm{D}}}\frac{1}{1/R_1^2+1/R_2^2}.$$

$${\tau }_{RC}^{\mathrm{s}}=\frac{{\lambda }_{\mathrm{D}}{R}_2}{D}\frac{1-R_1/R_2}{R_1/R_2+R_2/R_1}.$$

$$C^{\mathrm{c}}=\frac{2\pi \epsilon \ell }{{\lambda }_{\mathrm{D}}}\frac{1}{1/R_1+1/R_2}.$$

$${\tau }_{RC}^{\mathrm{c}}=\frac{{\lambda }_{\mathrm{D}}{R}_2}{D}\frac{\ln (R_2/R_1)}{1+R_2/R_1}.$$

$$\frac{2}{r^d}{\partial }_r[r^d{\partial }_r\varphi ]$$

$$\frac{{\partial }_{t}q}{D}=-\frac{1}{r^d}{\partial }_r[r^d{\partial }_{r}q]+{\kappa }^{2}q,$$

$$q(r,t=0)$$

$$\varphi (R_2,t>0)-\varphi (R_1,t>0)$$

$$J_q(R_1,t)=J_q(R_2,t)$$

$$q(r,t)$$

$$\frac{2}{r^d}{\partial }_r\left[r^d{\partial }_r\hat{\varphi }\right]$$

$$\frac{1}{r^d}{\partial }_r\left[r^d{\partial }_r\hat{q}\right]$$

$$\hat{\varphi }(R_2,s)-\hat{\varphi }(R_1,s)$$

$$\displaystyle-\partial_{r}\hat{q}-2\partial_{r}\hat{\phi}\Big{|}_{r=\{R_{1},R_{2}\}}$$

$$-2r^d{\partial }_r\hat{\varphi }=R_1^d{\partial }_r\hat{q}(R_1)+{\int }_{R_1}^r\mathrm{d}r{r}^d{\kappa }^2\hat{q},$$

$$-2r^d{\partial }_r\hat{\varphi }=R_2^d{\partial }_r\hat{q}(R_2)+{\int }_{R_2}^r\mathrm{d}r{r}^d{\kappa }^2\hat{q}.$$

$$0$$

$$-\frac{4\Delta \Phi }{s}$$

$$+\left[R_1^d{\partial }_r\hat{q}(R_1)+R_2^d{\partial }_r\hat{q}(R_2)\right]{\int }_{R_1}^{R_2}\frac{\mathrm{d}r}{r^d},$$

$${\hat{q}}^{\text{s}}$$

$${\Gamma }^{\text{s}}$$

$$\displaystyle\quad-\sinh[m(\xi-\bar{r})]+\sinh[m(1-\bar{r})]\big{\}}/\bar{r}\,,$$

$${{\rm Y}}^{\text{s}}$$

$$-\left(m^2-n^2-\frac{1}{\xi }\right)(1-\xi )\sinh [m(1-\xi )],$$

$${\hat{q}}^{\text{c}}\equiv \frac{2\Delta \Phi }{s}\frac{{\Gamma }^{\text{c}}}{{{\rm Y}}^{\text{c}}},$$

$$\displaystyle\Gamma^{\textrm{c}}=m^{2}\big{\{}\left[\xi K_{1}\!\left(m\xi\right)-K_{1}\!\left(m\right)\right]I_{0}\!\left(m\bar{r}\right)$$

$$\displaystyle\quad\quad\quad\quad+\left[\xi I_{1}\!\left(m\xi\right)-I_{1}\!\left(m\right)\right]K_{0}\!\left(m\bar{r}\right)\big{\}}\,,$$

$${{\rm Y}}^{\text{c}}=m(m^2-n^2)\xi \ln \xi \left[I_1\left(m\right)K_1\left(m\xi \right)-I_1\left(m\xi \right)K_1\left(m\right)\right]$$

$$\displaystyle\quad\quad+\frac{2n^{2}}{m}-n^{2}\big{[}\xi I_{1}\!\left(m\xi\right)K_{0}\!\left(m\right)+I_{1}\!\left(m\right)K_{0}\!\left(m\xi\right)$$

$$\displaystyle\quad\quad\quad\quad\quad\quad+\xi I_{0}\!\left(m\right)K_{1}\!\left(m\xi\right)+I_{0}\!\left(m\xi\right)K_{1}\!\left(m\right)\big{]}\,,$$

$${\Gamma }_{\text{eq}}^{\text{s}}$$

$$\displaystyle\quad-\sinh[n(\xi-\bar{r})]+\sinh[n(1-\bar{r})]\big{\}}/\bar{r}\,,$$

$${{\rm Y}}_{\text{eq}}^{\text{s}}$$

$$+2n,$$

$${\Gamma }_{\text{eq}}^{\text{c}}$$

$$+\left[\xi I_1\left(n\xi \right)-I_1\left(n\right)\right]K_0\left(n\bar{r}\right),$$

$${{\rm Y}}_{\text{eq}}^{\text{c}}$$

$$-\xi I_0\left(n\right)K_1\left(n\xi \right)-I_0\left(n\xi \right)K_1\left(n\right)+\frac{2}{n},$$

$$2\left(m-\frac{n^2}{m}\right)$$

$$+O\left\{\frac{n^2}{m}\exp [m(\xi -1)]\right\}.$$

$$q_1(r,t)$$

Full text

I Introduction

Many functionalities in nature and technology rely on the out-of-equilibrium behavior of electrolytes. Transport of ions through nerve membranes, for example, underlies the firing of neurons Hodgkin and Huxley (1952). Similarly, ionic fluxes in nanoporous carbon electrodes determine the power of supercapacitors Béguin et al. (2014) and the operation speed of capacitive deionization devices Suss et al. (2015). To optimize the performance of both capacitive devices through rational design, one needs a fundamental understanding of what sets the characteristic timescale $\tau$ of ionic response to electrode potentials. Since ions forming the electric double layer (EDL) must be partially drawn from a reservoir, $\tau$ could depend on “long” length scales like the electrode separation. This makes predicting $\tau$ with molecular simulations difficult, as typical simulation domains only capture a small portion of the nanoporous electrode structure Péan et al. (2014) or rely on simplified geometries Breitsprecher et al. (2018).

Theoretical predictions for $\tau$ typically concern either nontrivial electrode morphologies, treated approximately de Levie (1963); Biesheuvel and Bazant (2010), or concern the simplest of geometries, i.e., electrolytes between parallel planar blocking electrodes. In the latter case, the ionic charge density reacts to small suddenly imposed (DC) electrostatic potential differences on the timescale $\lambda_{\mathrm{D}}L/D$, which was derived with both microscopic and equivalent circuit model calculations Bazant et al. (2004); Janssen and Bier (2018). The bulk diffusion timescale $L^{2}/D$ can also appear, for instance, when large potentials are applied Bazant et al. (2004) or when the ionic diffusivities are unequal Alexe-Ionescu et al. (2007); Balu and Khair (2018). Different timescales even (with fractional powers of $\lambda_{\mathrm{D}}$ and $L$) appear when nonblocking electrodes are driven with an AC voltage Rubinstein et al. (2009). The analytical parallel-plate results of Refs. Bazant et al. (2004); Janssen and Bier (2018); Alexe-Ionescu et al. (2007); Balu and Khair (2018); Rubinstein et al. (2009) suggest that, to find $\tau$ for capacitive devices with complex nanoporous electrodes, one should identify relevant length scales and their relative importance. However, there are no general principles yet on how one should go about this task. Hence, it is timely to diminish the gap between analytical and molecular simulation predictions of electrolyte relaxation.

As a step towards an analytical understanding of the influence of nontrivial electrode morphology on ionic relaxation, in this article I discuss EDL capacitors with blocking concentric spherical or coaxial cylindrical electrodes (see Fig. 1). I use superscript s and c throughout this article to specify observables to either geometry. For both systems, the electrodes have radii $R_{1}$ and $R_{2}$ ($\Delta R=R_{2}-R_{1}>0$), respectively. The length $\ell$ of the cylinders is sufficiently large that I can ignore edge effects. While both systems then contain one relevant geometric length scale more than the parallel plate geometry (depending only on $L$) I will show that they allow for similar analytical insight. I assume spherical or axial symmetry in either case. Hence, all observables only depend on the radial distance $r$, with $R_{1}\leqslant r\leqslant R_{2}$. In between the electrodes is a dilute 1:1 electrolyte of dielectric constant $\varepsilon$. The ionic charge density, the difference between cationic and anionic densities, vanishes throughout the cell initially. Application of a small dimensionless potential difference $\Delta\Phi\ll 1$ (with electrostatic potentials measured in units of the thermal voltage $k_{\rm B}T/e$, with $k_{\rm B}T$ being the thermal energy and $e$ the proton charge) then drives the formation of EDLs at both electrode surfaces. Their equilibrium width is set by the Debye length $\lambda_{\mathrm{D}}=\kappa^{-1}$.

II $RC$ reasoning

Equivalent circuit representations of both setups in Fig. 1 contain two capacitors representing the EDLs at both electrode surfaces, a resistor for the electrolytic resistance, and a voltage source, all connected in series. In the spirit of Helmholtz, I treat the EDLs as dielectric capacitors of width $\lambda_{\mathrm{D}}$ and permittivity $\varepsilon$. Then, using that the capacitance of a dielectric capacitor of two conducting concentric spheres at $r_{1}$ and $r_{2}$ is $C^{\rm s}=4\pi\varepsilon/(1/r_{1}-1/r_{2})$, I find the capacitance of the EDL at the inner electrode ($r_{1}=R_{1},r_{2}=R_{1}+\lambda_{\mathrm{D}}$) as $C^{\rm s}_{R_{1}}\approx 4\pi\varepsilon R_{1}^{2}/\lambda_{\mathrm{D}}$, where I assumed $\lambda_{\mathrm{D}}/R_{1}\ll 1$. Likewise, the EDL at the outer electrode ($r_{1}=R_{2}-\lambda_{\mathrm{D}},r_{2}=R_{2}$) has a capacitance $C^{\rm s}_{R_{2}}\approx 4\pi\varepsilon R_{2}^{2}/\lambda_{\mathrm{D}}$ if $\lambda_{\mathrm{D}}/R_{2}\ll 1$. The two in-series EDLs have a total capacitance

[TABLE]

The resistance of the electrolyte is $R^{\rm s}=\rho/[4\pi(1/R_{1}-1/R_{2})]$, with $\rho=\lambda_{\mathrm{D}}^{2}/(\varepsilon D)$ being its resistivity Bazant et al. (2004). Multiplying $R^{\rm s}=\lambda_{\mathrm{D}}^{2}\Delta R/(4\pi\varepsilon DR_{1}R_{2})$ by $C^{\rm s}$ then yields the $RC$ time

[TABLE]

For the cylindrical electrode system, starting from $C^{\rm c}=2\pi\epsilon\ell/\ln(r_{2}/r_{1})$ and applying the same steps gives

[TABLE]

With the resistance $R^{\rm c}=\rho\ln(R_{2}/R_{1})/(2\pi\ell)$, I now find

[TABLE]

When $R_{1}\to R_{2}$, the electrodes locally resemble parallel plates and the relaxation times reduce to the familiar $\tau^{\rm c}_{RC}\approx\tau^{\rm s}_{RC}\approx\lambda_{\mathrm{D}}\Delta R/(2D)$. Conversely, for $R_{1}\ll R_{2}$, $\tau^{\rm s}_{RC}\approx R_{1}\lambda_{\mathrm{D}}/D$: The relaxation then only depends on the shortest geometric length scale. However, for general cases, $\tau^{\rm s}_{RC}$ and $\tau^{\rm c}_{RC}$ explicitly depend on both $R_{2}$ and $R_{1}$.

III Microscopic model

III.1 Governing equations

The dimensionless electrostatic potential $\phi(r,t)$ is related to the dimensionless ionic charge density $q(r,t)$ via the Poisson equation

[TABLE]

where ${d}=0$ for rectangular, ${d}=1$ for cylindrical, and ${d}=2$ for spherical coordinates. Moreover, $q(r,t)$ satisfies a continuity equation $\partial_{t}q=-\nabla\cdot\textbf{J}_{q}$. At small potentials $\phi(r)\ll 1$, $\textbf{J}_{q}=\hat{e}_{r}J_{q}$ with $J_{q}=-D(\partial_{r}q+2\partial_{r}\phi)$ Bazant et al. (2004); Janssen and Bier (2018). Inserting $\textbf{J}_{q}$ into the continuity equation and using Eq. (5) yields the Debye-Falkenhagen equation Debye and Falkenhagen (1928)

[TABLE]

subject to

[TABLE]

which account for initial charge neutrality, the suddenly imposed potential difference, and the no-flux (blocking) boundary conditions.

III.2 Solution to Laplace-transformed Debye-Falkenhagen equation

I determine $q(r,t)$ as follows. With Laplace transformations, the partial differential equation for $q(r,t)$ [Eq. (6)] turns into a solvable ordinary differential equation [Eq. (9b)] for its Laplace transformed counterpart $\hat{q}(r,s)=\mathcal{L}\{q(r,t)\}$ [likewise $\hat{\phi}(r,s)=\mathcal{L}\{\phi(r,t)\}$]. Then $q(r,t)$ is determined through

[TABLE]

with $s_{j}$ being the poles of $\hat{q}(r,s)$, labeled with $j$.

Applying $\mathcal{L}\{\,\}$ on both sides of Eqs. (5)–(7), I find

[TABLE]

with $k^{2}=\kappa^{2}+s/D$, subject to

[TABLE]

where I used Eq. (7a) for Eq. (9b). The solution to Eq. (9b) reads $\hat{q}^{\textrm{s}}(r)=a_{2}\exp{[-kr]}/r+b_{2}\exp{[kr]}/r$ for ${d}=2$ and $\hat{q}^{\textrm{c}}(r)=a_{1}I_{0}(kr)+b_{1}K_{0}(kr)$ for ${d}=1$, with $I_{0}\!$ and $K_{0}\!$ being modified Bessel functions of the first and second kind, respectively. The constants $a_{1}$, $a_{2}$, $b_{1}$, and $b_{2}$ could be fixed with the boundary conditions Eq. (10), which, however, inconveniently contain both $\hat{q}$ and $\hat{\phi}$. Aiming at two constraints on $\hat{q}$ only, I integrate Eq. (9a) over $\int_{R_{1}}^{r}\mathrm{d}r$ and use Eq. (10b) to find

[TABLE]

with $\partial_{r}\hat{q}(R_{1})$ shorthand for \partial_{r}\hat{q}(r)\big{|}_{r=R_{1}}.

Repeating the same calculation for $\int_{R_{2}}^{r}\mathrm{d}r$ gives

[TABLE]

The difference and the sum (integrated over $\int_{R_{1}}^{R_{2}}\mathrm{d}r$) of Eqs. (11) and (12) read

[TABLE]

which are two constraints on $\hat{q}^{\textrm{s}}(r,s)$ and $\hat{q}^{\textrm{c}}(r,s)$ each, which fix the constants $a_{1}$, $a_{2}$, $b_{1}$, and $b_{2}$ therein 111For ${d}=0$, inserting $\hat{q}(r)=a_{0}\exp{[-kr]}+b_{0}\exp{[kr]}$ gives $\hat{q}(r)$ as reported in Refs. Bazant et al. (2004); Janssen and Bier (2018).. For spherical electrodes, I find

[TABLE]

while for cylindrical electrodes I find

[TABLE]

where $m\equiv kR_{2}$, $n\equiv\kappa R_{2}$, $\xi\equiv R_{1}/R_{2}$, and $\bar{r}=r/R_{2}$. Here $n$ measures the thickness of the EDLs relative to the system size. For most practical devices, the nanometer-sized EDLs are well separated, i.e., $n\gg 1$.

III.3 Equilibrium

The pole $s_{0}\equiv 0$ in Eqs. (14a) and (15a) sets the equilibrium charge density through $q_{\textrm{eq}}(r)\equiv\text{Res}\left(\hat{q}(r,s),0\right)$. This amounts to $q^{\textrm{s}}_{\textrm{eq}}(r)\equiv 2\Delta\Phi\Gamma_{\textrm{eq}}^{\textrm{s}}/\Upsilon_{\textrm{eq}}^{\textrm{s}}$, with

[TABLE]

and $q^{\textrm{c}}_{\textrm{eq}}(r)\equiv 2\Delta\Phi\Gamma^{\textrm{c}}_{\textrm{eq}}/\Upsilon^{\textrm{c}}_{\textrm{eq}}$, with

[TABLE]

respectively [i.e., $\Gamma_{\textrm{eq}}\equiv\Gamma(m=n)$ and $\Upsilon_{\textrm{eq}}\equiv\Upsilon(m=n)$]. Note that $q^{\textrm{s}}_{\textrm{eq}}(r)$ and $q^{\textrm{c}}_{\textrm{eq}}(r)$ can be derived more easily: At equilibrium, $J_{q}(r)=0$ gives $q(r)=-2\phi(r)$, and thus $\partial_{r}[r^{2}\partial_{r}\phi]=r^{2}\kappa^{2}\phi$ [Eq. (5)]. The solution to this equation, $\phi^{\textrm{s}}(r)=A_{2}\exp{[-\kappa r]}/r+B_{2}\exp{[\kappa r]}/r$, contains two constants ($A_{2}$ and $B_{2}$), which are fixed with Eq. (7b) and particle conservation, \int_{R_{1}}^{R_{2}}\mathrm{d}r\,r^{2}q=0\Rightarrow r^{2}\partial_{r}\phi\big{|}_{R_{1}}^{R_{2}}=0. With $q(r)=-2\phi(r)$, $q^{\textrm{s}}_{\textrm{eq}}(r)$ then trivially follows. The same steps for $d=1$ lead to $q^{\textrm{c}}_{\textrm{eq}}(r)$. Further details on the equilibrium EDL near curved electrodes can be found in Refs. Lian et al. (2016); Reindl et al. (2017a); Bossa et al. (2019); Asta et al. (2019) and in the Appendix, where I confirm Eqs. (1) and (3) with $q^{\textrm{s}}_{\textrm{eq}}$ and $q^{\textrm{c}}_{\textrm{eq}}$, respectively.

III.4 Relaxation time

For the relaxation of $q(r,t)$, I need to determine the locations of the poles $s_{j}\in\mathbb{C}$ of $\hat{q}(r,s)$. However, instead of immediately focusing on $s_{j}$, Ref. Janssen and Bier (2018) showed for the case of planar electrodes that it is easier to determine the corresponding poles $m_{j}\in\mathbb{C}$ of $\hat{q}(r,m)$ first. As in Ref. Janssen and Bier (2018), each $m_{j}$ that I could find was either purely real or purely imaginary (discussed below). Thus, all corresponding $s_{j}=(m_{j}^{2}-n^{2})D/R_{2}^{2}$ are real and, as is turns out, $s_{j\neq 0}<0$. As I am interested in the late-time response of $q(r,t)$, I focus here on $s_{1}$, the pole closest to $s_{0}$, as this pole sets the late-time relaxation time $\tau_{1}=-1/s_{1}$.

To find $s_{1}$, I first note that neither $\Gamma^{\textrm{s}}(m)$ nor $\Gamma^{\textrm{c}}(m)$ has poles in $m\in\mathbb{C}$. Thus, all $s_{j\neq 0}$ come from the zeros of $\Upsilon^{\textrm{s}}(m)$ and $\Upsilon^{\textrm{c}}(m)$, respectively. Both $\Upsilon^{\textrm{s}}(m)$ and $\Upsilon^{\textrm{c}}(m)$ oscillate around zero on the imaginary $m$-axis [Fig. 2(a)] and hence contribute to Eq. (8) with infinitely many poles. However, only the zero at the smallest $m$-value has the potential of leading to $s_{1}$; all zeros further along the imaginary $m$-axis give smaller $s_{j}$, hence faster decaying modes.

Figures 2(b) and 2(c) show $\Upsilon(m)$ for $m\in\mathbb{R}$ at $n=15$ and $n=100$, respectively. For certain values of $\xi$, solutions $m_{1}=\mathfrak{M}$ to $\Upsilon^{\textrm{s}}(m)=0$ and $\Upsilon^{\textrm{c}}(m)=0$ are visible there. For $n=15$, $\mathfrak{M}$ disappears at $\xi=0.767$ ($\Upsilon^{\textrm{s}}$) and $\xi=0.768$ ($\Upsilon^{\textrm{c}}$), respectively. Conversely, at $n=100$, $\mathfrak{M}$ persists until higher $\xi$, disappearing only at $\xi=0.965$ (both $\Upsilon^{\textrm{s}}$ and $\Upsilon^{\textrm{c}}$). In both cases, $m_{1}$ disappears through the origin (not shown), after which it reappears as a zero $m_{1}=i\mathcal{M}$ (with $\mathcal{M}\in\mathbb{R}$) on the imaginary $m$-axis [see the inset of Fig. 2(b)] that moves away from the origin with increasing $\xi$ 222A similar transition from a solution $\mathfrak{M}$ to $\mathcal{M}$ occurs for planar electrodes Janssen and Bier (2018). In that simpler geometry (without $\xi$), the transition happens at $n=\sqrt{3}$ always. (This transition from $\mathfrak{M}$ to $\mathcal{M}$ occurs also at small $\xi$ and $n$ with decreasing $\xi$.) Associated with these zeros, $s_{1}$ is either $s_{1}=\left(\mathfrak{M}^{2}-n^{2}\right)D/R_{2}^{2}$ or $s_{1}=-\left(\mathcal{M}^{2}+n^{2}\right)D/R_{2}^{2}$.

Figure 3 shows $\tau_{1}=-1/s_{1}$ for both setups at $n=\kappa R_{2}=15$ (dashed lines) and $n=100$ (solid lines). At red crosses, $m_{1}$ transitions from $\mathfrak{M}$ to $i\mathcal{M}$. The opposite transition (from $\mathcal{M}$ to $\mathfrak{M}$), indicated with orange boxed crosses, does not occur at $n=100$ and occurs for $\Upsilon^{\textrm{c}}$ only around $\xi=10^{-12}$. The plateaus at small $\xi$ are understood as follows. For $\xi\ll 1$, $\Upsilon^{\textrm{s}}=0$ reduces to $\tanh m=m$, whose only solution $m_{1}=0$ gives $\tau D/(R_{2}\lambda_{\mathrm{D}})\approx 1/n$. Also shown in Fig. 3 with circles are the equivalent circuit model predictions $\tau^{\rm s}_{RC}$ [Eq. (2)] and $\tau^{\rm c}_{RC}$ [Eq. (4)]. For $\xi>0.1$ and $n=15$, $\tau_{1}$ and $\tau_{RC}$ are qualitatively similar. For larger $n$, $\tau_{1}$ and $\tau_{RC}$ become identical. This is understood analytically as follows. For $m\in\mathbb{R}$ and $n\gg 1$, $\Upsilon^{\textrm{s}}(m)=0$ reduces to

[TABLE]

From Fig. 2(c) I see that $\mathfrak{M}\approx n$ if $n\gg 1$. Inserting the approximation $\mathfrak{M}_{\textrm{ap}}=n-\delta$ into Eq. (III.4) and keeping terms up to $O(\delta)$, I find $\delta=(\xi+1/\xi)/[2(1-\xi)]$. This reproduces Eq. (2): $\tau^{\textrm{s}}_{1}\approx R_{2}\lambda_{\mathrm{D}}(1-\xi)/[D(\xi+1/\xi)]$.

Similarly, for $n\gg 1$, $\Upsilon^{\textrm{c}}(m)=0$ amounts to $m\left(m^{2}-n^{2}\right)\xi\ln{\xi}=n^{2}\left(1+\xi\right)+O\left\{(\exp[2n(\xi-1)]\right\}$. [I used Hankel’s large argument expansions $I_{\alpha}(z)\sim\exp{(z)}/\sqrt{2\pi z}$ and $K_{\alpha}(z)\sim\sqrt{\pi}\exp{(-z)}/\sqrt{2z}$ here, which imply $I_{\alpha}(n)\gg I_{\alpha}(n\xi)$, $K_{\alpha}(\xi n)\gg K_{\alpha}(n)$, $I_{1}(n)/I_{0}(n)\to 1$, and $K_{1}(n)/K_{0}(n)\to 1$]. Inserting $\mathfrak{M}_{\textrm{ap}}=n-\epsilon$ and keeping terms up to $O(\epsilon)$ yields $\epsilon=-(1+\xi)/(2\xi\ln{\xi})$. This gives $\tau^{\textrm{c}}_{1}\approx-R_{2}\lambda_{\mathrm{D}}\xi\ln\xi/[D(1+\xi)]$, i.e., Eq. (4).

III.5 Relaxation of $q(r,t)$

I use that, close to $s_{1}$, $\Upsilon(s)\overset{s\to s_{1}}{=}\partial_{m}\Upsilon(m_{1})\times(s-s_{1})R_{2}^{2}/(2mD)$. The slowest relaxation mode $q_{1}(r,t)\equiv\text{Res}\left(\hat{q}(r,s)\exp(st),s_{1}\right)$ now amounts to

[TABLE]

Explicit expressions for $\partial_{m}\Upsilon(m)$ for the respective geometries follow with Eqs. (14c) and (15c). I find

[TABLE]

and

[TABLE]

Truncating the sum in Eq. (8) after $j=1$, I approximate the relaxation of $q(r,t)$ by

[TABLE]

In Fig. 4, I compare $q_{\textrm{ap}}(r,t)$ to numerical inversions of Eqs. (14) and (15) with the ’t Hoog algorithm, respectively. I observe a stronger asymmetry in $q^{\textrm{s}}$ than in $q^{\textrm{c}}$, which must stem from the difference in their ratios of inner to outer electrode surface areas ($\xi^{2}$ and $\xi$ in either case, respectively). Note that, at $n=15$ as used here, the numerically determined $\mathfrak{M}^{\textrm{s}}=12.15$ and $\mathfrak{M}^{\textrm{c}}=12.67$ deviate substantially from their analytical approximations $\mathfrak{M}^{\textrm{s}}_{\textrm{ap}}=12.58$ and $\mathfrak{M}^{\textrm{c}}_{\textrm{ap}}=13.10$, respectively; hence, I use the former. Clearly, $q^{\textrm{s}}_{\textrm{ap}}(r,t)$ describes $q^{\textrm{s}}(r,t)$ well for $t\approx\tau^{\textrm{s}}_{1}=0.0129R^{2}_{2}/D$, while deviations are visible at early times. The same is true for $q^{\textrm{c}}_{\textrm{ap}}(r,t)$ and $q^{\textrm{c}}(r,t)$ for times around $\tau^{\textrm{c}}_{1}D/R^{2}_{2}=0.0155$. Better capturing $q(r,t)$ at early times requires truncating the sum in Eq. (8) at higher $j$, which I leave for future work.

Finally, as I focused on electrolytes with equal ionic diffusivities $D_{+}=D_{-}$ in this article, the results derived here can be expected to be accurate for, for instance, KCl, RbBr, and CsBr (which have $D_{+}/D_{-}=0.97,1.00$, and $1.00$, respectively Agar et al. (1989)), but less so for NaCl ($D_{+}/D_{-}=0.66$). Electrolytes with different ionic diffusivities will probably relax on two timescales: a fast $RC$ timescale as described here and a slower diffusive timescale that becomes more important the more $D_{+}$ and $D_{-}$ differ Alexe-Ionescu et al. (2007); Balu and Khair (2018); Janssen and Bier (2019). Finding the precise functional form of this diffusive timescale can be done with calculations along the lines of the ones presented in Refs. Balu and Khair (2018); Janssen and Bier (2019).

IV Conclusions

I have studied the influence of electrode morphology on the relaxation of EDL capacitors, both with equivalent circuit models and with the (microscopic) Debye-Falkenhagen equation. For two different curved-electrode geometries, I have shown that the timescale of ionic response to an applied electrostatic potential explicitly depends on all geometric length scales present. The uplifting message is that, for thin EDLs (a case of high practical relevance), easily obtainable $RC$ times capture the ionic relaxation times decently. Conversely, for thick EDLs, corrections must be taken into account. These results form a small step towards an analytical understanding of the relaxation of supercapacitors and deionization devices; however, I expect complications at each further step of the way from planar to nanoporous electrodes.

Acknowledgements.

I gratefully acknowledge a discussion with Aymar de Lichtervelde and Pedro de Souza on Ref. de Lichtervelde et al. (2019), which inspired me to write this article. I thank Martin Z. Bazant and Watse Sybesma for insightful comments and I acknowledge S. Dietrich for support.

Appendix A Equilibrium surface charge and capacitance

The total surface charge $Q^{\textrm{s}}_{1}$ of the smaller electrode is related to its unit surface charge density $\sigma^{\textrm{s}}_{1}$ by $Q^{\textrm{s}}_{1}=4\pi R_{1}^{2}e\sigma^{\textrm{s}}_{1}$. With Gauss’s law $e^{2}\sigma^{\textrm{s}}_{1}=-\varepsilon k_{\rm B}T\partial_{r}\phi(R_{1})$ and $\phi(r)=-q(r)/2$ I find $Q_{1}e=2\pi R_{1}^{2}\varepsilon k_{\rm B}T\partial_{r}q^{\textrm{s}}_{1}(R_{1})$. Likewise, the charge on the larger electrode reads $Q^{\textrm{s}}_{2}e=-2\pi R_{2}^{2}\varepsilon k_{\rm B}T\partial_{r}q^{\textrm{s}}(R_{2})$. With Eq. (16), it is now straightforward to verify that the electrodes carry opposite surface charge $Q^{\textrm{s}}_{1}=-Q^{\textrm{s}}_{2}$.

Likewise, the capacitance $C=eQ_{2}/(k_{\rm B}T\Delta\Phi)$ reads

[TABLE]

For the cylindrical electrode system, the surface charge $Q^{\textrm{c}}_{1}=2\pi R_{1}\ell e\sigma^{\textrm{c}}_{1}$ of the inner electrode amounts to $Q^{\textrm{c}}_{1}e=\pi R_{1}\ell\varepsilon k_{\rm B}T\partial_{r}q(R_{1})$. Similarly, $Q^{\textrm{c}}_{2}e=-\pi R_{2}\ell\varepsilon k_{\rm B}T\partial_{r}q(R_{2})$. Again, with Eq. (17), $Q^{\textrm{c}}_{1}=-Q^{\textrm{c}}_{2}$ can be shown to hold. The capacitance reads

[TABLE]

where, going to the third line, I again used Hankel’s large argument expansion, stated below Eq. (III.4). Note that Eqs. (A) and (A) are equivalent to Eqs. (1) and (3), respectively.