Current Distribution on Capacitive Electrode-Electrolyte Interfaces

TL;DR

This paper presents an analytical model describing how current distribution evolves on capacitive electrode-electrolyte interfaces over time, with experimental validation, relevant for electrochemical and neural interface design.

Contribution

It introduces a new analytical model for current distribution dynamics at capacitive interfaces, validated experimentally, enhancing understanding of electrochemical interface behavior.

Findings

Steady current distribution is proportional to capacitance per unit area.

Non-uniform initial current distribution becomes uniform over time.

Model applies to electrode geometries like disks with experimental support.

Abstract

The distribution of electric current on an electrode surface in electrolyte varies with time due to charge accumulation at a capacitive interface, as well as due to electrode kinetics and concentration polarization in the medium. Initially, the potential at the electrode-electrolyte interface is uniform, resulting in a non-uniform current distribution due to the uneven ohmic drop of the potential in the medium. Over time, however, the non-uniform current density causes spatially varying rate of the charge accumulation at the interface, breaking down its equipotentiality. We developed an analytical model to describe such transition at a capacitive interface when the current is below the mass-transfer limitation, and demonstrated that the steady distribution of the current is achieved when the current density is proportional to the capacitance per unit area, which leads to linear voltage ramp at the electrode. More specific results regarding the dynamics of this transition are provided for a disk electrode, along with an experimental validation of the theoretical result. These findings are important for many electrochemical applications, and in particular, for proper design of the electro-neural interfaces.

Equations

$$\phi (\mathbf{r},t):=\Phi (\mathbf{r},t),\mathbf{r}\in A.$$

$$\phi (\mathbf{r},t)+U(\mathbf{r},t)=V(t).$$

$$\mathbf{i}(\mathbf{r},t):=-\frac{1}{\rho }\nabla \Phi (\mathbf{r},t),$$

$$\mathbf{i}(\mathbf{r})\cdot \mathbf{n}(\mathbf{r})=0,\mathbf{r}\in D,$$

$$i(\mathbf{r},t):=\mathbf{i}(\mathbf{r},t)\cdot \mathbf{n}(\mathbf{r}),\mathbf{r}\in A.$$

$$i_a(\mathbf{r},t)+i_b(\mathbf{r},t)=i(\mathbf{r},t).$$

$$i_a(\mathbf{r},t)=C_p(\mathbf{r}){\dot{U}}_p(\mathbf{r},t),$$

$$i_b(\mathbf{r},t)=C_d(\mathbf{r})\dot{U}(\mathbf{r},t).$$

$$R_{ct}(\mathbf{r})i_a(\mathbf{r},t)+U_p(\mathbf{r},t)=U(\mathbf{r},t).$$

$$\Delta \Phi =0,$$

$$\hat{\mathbf{S}}\phi =i.$$

$$i_a=\hat{\mathbf{S}}\phi -C_d(V^{\prime }-\dot{\phi }).$$

$$R_{ct}{\dot{i}}_a+C_p^{-1}{i}_a=V^{\prime }-\dot{\phi }.$$

$$\begin{array}{rl}& {\hat{\mathbf{R}}}_a{R}_{ct}{C}_p{C}_d\ddot{\phi }+\left[{\hat{\mathbf{R}}}_a{C}_d+\left({\hat{\mathbf{R}}}_a+R_{ct}\right)C_p\right]\dot{\phi }+\phi \\ =& {\hat{\mathbf{R}}}_a{R}_{ct}{C}_p{C}_d{V}^{\prime \prime }+{\hat{\mathbf{R}}}_a(C_d+C_p)V^{\prime }.\end{array}$$

$$\ddot{\phi }=\dot{\phi }=0.$$

$$i(\mathbf{r})=\hat{\mathbf{S}}\phi =R_{ct}{C}_p{C}_d{V}^{\prime \prime }(t)+(C_d+C_p)V^{\prime }(t).$$

$$R_{ct}\frac{C_p{C}_d}{C_d+C_p}{V}^{\prime \prime }+V^{\prime }=v.$$

$$i(\mathbf{r})=\left[C_d(\mathbf{r})+C_p(\mathbf{r})\right]v.$$

$$p(x)=x^2{\hat{\mathbf{R}}}_a{R}_{ct}{C}_p{C}_d+x\left[{\hat{\mathbf{R}}}_a{C}_d+\left({\hat{\mathbf{R}}}_a+R_{ct}\right)C_p\right]+\hat{I},$$

$$\hat{\mathbf{P}}:=-{\hat{\mathbf{R}}}_a{\left({\hat{\mathbf{R}}}_a+R_{ct}\right)}^{-1},$$

$$p(x)\approx \tilde{p}(x):=p(x)+x{\hat{\mathbf{R}}}_a{C}_d\hat{\mathbf{P}}.$$

$$\begin{array}{r}\tilde{p}(x)=\left(x{\hat{\mathbf{T}}}_b+\hat{I}\right)\left(x{\hat{\mathbf{T}}}_a+\hat{I}\right).\end{array}$$

$$\begin{array}{rl}& \left({\hat{\mathbf{T}}}_b{\partial }_t+\hat{I}\right)\left({\hat{\mathbf{T}}}_a{\partial }_t+\hat{I}\right)\phi \\ =& \left[\left({\hat{\mathbf{T}}}_b{\partial }_t+\hat{I}\right)\left({\hat{\mathbf{T}}}_a{\partial }_t+\hat{I}\right)-R_{ct}{C}_p{\partial }_t-\hat{I}\right]V.\end{array}$$

$$\begin{array}{rl}\phi =& {\left(j\omega {\hat{\mathbf{T}}}_a+\hat{I}\right)}^{-1}{\left(j\omega {\hat{\mathbf{T}}}_b+\hat{I}\right)}^{-1}\\ & \left[-{\omega }^2{\hat{\mathbf{T}}}_b{\hat{\mathbf{T}}}_a+j\omega {\hat{\mathbf{R}}}_a(C_p+C_d)\right]V_0{e}^{j\omega t}.\end{array}$$

$$\begin{array}{rl}\hat{\mathbf{Z}}=V{i}^{-1}=& {\left[-{\omega }^2{\hat{\mathbf{T}}}_b{\hat{\mathbf{T}}}_a+j\omega {\hat{\mathbf{R}}}_a(C_p+C_d)\right]}^{-1}\\ & \left(j\omega {\hat{\mathbf{T}}}_b+\hat{I}\right)\left(j\omega {\hat{\mathbf{T}}}_a+\hat{I}\right){\hat{\mathbf{R}}}_a,\end{array}$$

$$\phi$$

$$i$$

$$\phi$$

$$i$$

$$\hat{\mathbf{Z}}\approx C_p^{-1}\hat{\mathbf{S}}\left({\hat{\mathbf{R}}}_a+R_{ct}\right)C_p{\hat{\mathbf{R}}}_a.$$

$${\hat{\mathbf{T}}}_a\dot{\phi }+\phi ={\hat{\mathbf{R}}}_a{C}_p\dot{V}$$

Full text

††preprint: APS/123-QED

I Introduction

Dynamics of the charge transfer across the electrode-electrolyte interfaces is of great importance in electrochemistry in general, and for many applications, including batteries, electroplating, electrolysis, chemical sensors and, in particular, bioelectronics. The distribution of current and voltage across such an interface is governed by multiple mechanisms, including the concentration polarization of the reactants in the medium, the kinetics of the electrode reactions, the ohmic drop in the bulk of electrolyte and charging of the electric double layer. The effects of the concentration polarization are modeled by the Warburg impedance, which is only significant at high current density when the reactant concentration is considerably affected by the mass-transfer limitationBard and Faulkner (2001); Newman and Thomas-Alyea (2012). The electrode kinetics is associated with faradaic electrochemical reactions, and modeled by the charge transfer resistance varying with voltage, which is considered linear at lower current density according to the Butler-Volmer modelBard and Faulkner (2001); Newman and Thomas-Alyea (2012). The access resistance – the ohmic drop in the medium – is determined only by the electrode geometry and electrolyte conductivity. The electric double layer is modeled as a capacitance, where the Helmholtz plane in the electrolyte serves as the ”plate” on the electrolyte side of this capacitive interfaceMyland and Oldham (2005). The double-layer capacitance is typically on the order of $10-$20\text{,}\mathrm{\SIUnitSymbolMicro F}\text{/}{\mathrm{cm}}^{2}$$ for inert materials including carbonBleda-Martínez et al. (2005), platinumPajkossy and Kolb (2001) and goldGore et al. (2010). Additionally, some materials can exhibit a range of quasi-continuous oxidation states, enabling reversible storage of much larger amount of charge than in a typical double-layer capacitance, and therefore known as pseudocapacitanceSugimoto et al. (2004); Soon and Loh (2007); Grupioni et al. (2002). Together, the double-layer capacitance and the pseudocapacitance are often called supercapacitanceHalper and Ellenbogen (2006).

NewmanNewman (1966a) calculated the primary current distribution at the interface of an equipotential (EP) disk electrode. It has been pointed out that both the electrode kinetics and the mass-transfer limitation result in secondary current distributions being more uniform than the primary oneNewman (1966b); West and Newman (1989). When the ohmic drop in the medium is the dominant part of the impedance, the calculation of the access resistance based on the EP boundary condition from Newman (1966a) is broadly usedAlbery et al. (1971); Weaver and Chizmadzhev (1996); Lasia (2002); Boinagrov et al. (2015). The transient charge redistribution within the double layer on a disk electrode made of the same material has been described for controlled potentialNisanciogğlu and Newman (1973a) and for controlled currentNisanciogğlu and Newman (1973b), respectively. The frequency dispersion of such interfaces was studied in Newman (1970).

However, there was no theory describing the current distribution on electrodes of all geometries or with multiple surface materials. Such theory is of interest in many electrochemical applications involving capacitive coupling electrodes, and especially for neural stimulation, where various geometries and materials are used in different applications. For neural stimulation electrodes, the distribution of the electric current affects the stimulation thresholds and tissue safety, and extra care should be taken to avoid irreversible electrochemical reactions. Therefore, materials of large charge storage capacity are often used to minimize the voltage swing, and to ensure that the charge transfer is fully reversible, i.e. pseudo-capacitive.

In this study we demonstrate, for any electrode geometry, that in the absence of the concentration polarization, the steady state current distribution is achieved when the current density is proportional to the surface capacitance per unit area (PCD), where the boundary condition is not necessarily EP and the electrode potential converges to a linear ramp. For an electrode made of the same material, PCD implies uniform current density (UCD). Initially, current begins to flow at a non-uniform density from the EP surface, but over time the uneven charge accumulation at the Helmholtz plane, as well as that in the pseudocapacitance, begins to affect the voltage drop across the double layer. Such uneven potential at the Helmholtz plane rearranges the electric field in the electrolyte, and hence redistributes the current density, until the system reaches the PCD steady state.

Note that we model the possibly changing equilibrium potential of the associated electrochemical reactions by the pseudocapacitance, while Nisanciogğlu and Newman (1973a, b) assumed constant reaction potentials. Therefore, in their notion, the steady state was resistive and unrelated to the surface capacitance. More recent studies Myland and Oldham (2005); Behrend et al. (2008) assumed ideally polarizable disk electrodes – with no faradaic reactions – and described the transition from the primary current distribution to the steady state with finite-element models, yielding only numerical solutions. However, without analytical description, fitting the numerical results to the RC approximations provides only a limited understanding of the transition, let alone that such finite-element models are usually intractable for an arbitrary electrode geometry. Nevertheless, the results have been widely adopted in practical applicationsCantrell et al. (2007); Myland and Oldham (2014); Bieniasz (2015); Sue et al. (2015).

The general model in this paper presents the transition from the initial current distribution to PCD for any electrode geometry and material composition, while considering the effects of the supercapacitance, the electrode kinetics and the ohmic drop. We develop a framework to study capacitive interfaces with sinusoidal waveforms, chronoamperometry (controlled potential) or chronopotentiometry (controlled current), with the bra-ket notation. We demonstrate the application of this framework to a disk electrode, which agrees with the previous analytical results. We also demonstrate validation of some of these results experimentally.

II The System Model

Typically, an equivalent circuit model of the electrode-electrolyte interface with supercapacitance includes the double-layer capacitance $C_{d}$, the pseudocapacitance $C_{p}$, the Faradaic leakage resistance $R_{f}$, the charge transfer resistance $R_{ct}$, the Warburg impedance $Z_{w}$, and the access resistance $R_{a}$Bard and Faulkner (2001); Conway and Pell (2003), as shown in FIG. 1(a). The complex Warburg impedance has a constant angle of $-\pi/4$, which is a result of the phase delay between the current and the concentration polarization, stemming from the diffusion of reactants. It is a special case of constant phase element (CPE). In more complicated models, more than one CPEs may be included, whose detailed mechanisms often remain unclearBarsoukov and Macdonald (2005).

We are interested in the conditions where the current density is below the mass-transfer limitation, and hence $Z_{w}$ is insignificant. The limiting current density is given by equation (28) of Newman (1966b) for a disk electrode. We also assume negligible Faradaic leakage. The charge transfer across the interface in this case is governed by the distributed capacitance and charge transfer resistance at the interface, and by the ohmic drop in the medium. Thereby, our circuit can be simplified to that shown in FIG. 1(b). The model studied in Myland and Oldham (2005) is a special case when $C_{p}=0$ and the one in Nisanciogğlu and Newman (1973a, b) is when $C_{p}\xrightarrow{}+\infty$.

For an extended electrode, the capacitance, the charge transfer resistance and the access resistance are distributed, as illustrated in FIG. 2. Each location on the interface is approximated by a discrete circuit in FIG. 1(b). Although illustrated with discrete components, our mathematical treatment makes no explicit discretization. Note that we only study half of the electrochemical cell, assuming that the current is collected on a large counter electrode infinitely far away.

In FIG. 2, $E$ denotes the subset of the 3-dimensional space occupied by the electrolyte. Its boundary with the electrode is denoted by $A$, and with an insulating surface, denoted by $D$. $A\bigcup D=\partial E\subset E$. Let $\Phi(\bm{r},t)$ denote the potential distribution in $E$ as a function of both the spatial variable $\bm{r}$ and time $t$. We choose $\Phi(\infty)=0$, and define $\varphi$ as the 2-dimensional restriction of $\Phi$ on $A$:

[TABLE]

The electrode is EP in its bulk, whose potential as a function of time is denoted as $V(t)$. A non-uniform potential drop across the surface, which also varies with time, is denoted by $U(\bm{r},t)$:

[TABLE]

$U$ is the voltage across $C_{d}$, and we denote the potential drop across $C_{p}$ by $U_{p}(\bm{r},t)$. Let $C_{p}(\bm{r})$ and $C_{d}(\bm{r})$ be the pseudocapacitance and the double-layer capacitance per unit area on $A$, respectively. $R_{ct}{(\bm{r})}$ is the charge transfer resistance times unit area on $A$. Typically,the pseudocapacitance is much larger than the double-layer capacitanceHalper and Ellenbogen (2006). For an area where no electrochemical reactions take place, we may set $C_{d}=0$, $R_{ct}=0$ and $C_{p}$ the double-layer capacitance, so that $C_{p}(\bm{r})\gg C_{d}(\bm{r})$ for $\forall\bm{r}\in A$.

In $E$, the current density is

[TABLE]

where $\rho$ is resistivity of the electrolyte. On $D$, insulation implies zero normal current:

[TABLE]

where the unitary $\bm{n}(\bm{r})$ is normal to the surface at $\bm{r}$ pointing to the electrolyte side. We choose the direction of $\bm{n}(\bm{r})$ as positive for current flow, and define the normal current density on $A$:

[TABLE]

Only the normal current component contributes to the charge accumulation on $A$. Conceptually, we may divide $i$ into two components: the faradaic component $i_{a}$ (the $C_{p}$–$R_{ct}$ path) and $i_{b}$ that charges the double layer (the $C_{d}$ path). By definition:

[TABLE]

$i_{a}$ drives the voltage of the pseudocapacitance by

[TABLE]

and similarly, for $i_{b}$ we have:

[TABLE]

The two current paths have the same potential drop:

[TABLE]

The potential $\Phi$ satisfies the Laplace’s equation:

[TABLE]

with the boundary conditions (1), (4) and $\Phi(\infty)\to 0$ as $\hbox{dist}(\bm{r},A)\to+\infty$. Given the boundary condition value $\varphi$, the potential $\Phi$ is fully determined, thus so is $i$. As (10) and all its boundary conditions are linear, there is a linear mapping from $\varphi$ to $i$, which we denote as a linear operator $\hat{\bm{S}}$:

[TABLE]

Because of the uniqueness of the electric field, $\hat{\bm{S}}$ is reversible. Define $\hat{\bm{R}}_{a}:=\hat{\bm{S}}^{-1}$. Combining (2), (6), (8) and (11), so that (6) gives

[TABLE]

We take the derivative of (9), and use (2) and (7):

[TABLE]

Combine (12) and (13):

[TABLE]

At the steady state of current distribution, by (11), $\varphi$ no longer changes with time:

[TABLE]

Thus at steady state we have

[TABLE]

In (16), the left hand side is time-independent, implying there is a constant $v$ so that:

[TABLE]

Therefore, $V^{\prime}$ converges to $v$ with time constant $R_{ct}{C_{p}C_{d}}/(C_{d}+C_{p})$, which may vary in space. The ”steady” steady state has $V^{\prime\prime}=0$, and yields

[TABLE]

Because $v$ is constant in space, at steady state the (normal) current density on the electrode surface is proportional to the capacitance per unit area. Note that a capacitive interface requires $v\neq 0$ for non-zero steady states of current distribution.

By (14), the characteristic polynomial of the system is

[TABLE]

where $\hat{I}$ is the identity operator. As proven in Appendix A, the operator $\hat{\bm{S}}$ is positive-definite, and so is $\hat{\bm{R}}_{a}$. $R_{ct}$ and $C_{d}$ are positive functions. Since $C_{p}\gg C_{d}$, the coefficient of the first-order term of (19) is dominated by $C_{p}$. Therefore, we may perturb the term $\hat{\bm{R}}_{a}C_{d}$ by its magnitude to make the system more tractable. We define a perturbation operator

[TABLE]

whose operator norm is bounded by 1. We perturb $p(x)$:

[TABLE]

Let $\hat{\bm{T}}_{a}=\left(\hat{\bm{R}}_{a}+R_{ct}\right)C_{p}$, and

$\hat{\bm{T}}_{b}=\hat{\bm{R}}_{a}R_{ct}C_{d}{\left(\hat{\bm{R}}_{a}+R_{ct}\right)}^{-1}$. We have:

[TABLE]

From (22), we know that the system has two sets of characteristic times, $\{\tau_{a}\}$ and $\{\tau_{b}\}$, corresponding to the eigenvalues of $\hat{\bm{T}}_{a}$ and $\hat{\bm{T}}_{b}$, respectively. Let $\{\Upsilon_{a}\}$ and $\{\Upsilon_{b}\}$ be the normalized (dimensionless) eigenfunctions of $\hat{\bm{T}}_{a}$ and $\hat{\bm{T}}_{b}^{\dagger}$, respectively. Each of the eigenfunctions corresponds to a eigenmode of $\varphi$, which is a potential distribution on $A$ that elicits a current such that $\varphi$ changes proportionally to itself. As mentioned, $\hat{\bm{R}}_{a}$ is positive-definite, so $\{\tau_{a}\}$ and $\{\tau_{b}\}$ are positive. Furthermore, if the surface is uniform, viz. $R_{ct}$, $C_{p}$ and $C_{d}$ are constant, $\hat{\bm{T}}_{a}$ and $\hat{\bm{T}}_{b}$ are positive-definite and thus the eigenmodes are orthogonal within each set. Because $C_{p}\gg C_{d}$, $\{\tau_{a}\}$ are much larger than $\{\tau_{b}\}$. Thus, $\max\{\tau_{a}\}$ is the dominant time constant, and $\{\varphi_{b}\}$ decay much faster than $\{\varphi_{a}\}$. Different eigenmodes of the same operator have different time constants because of the shape of the eigenmodes. Intuitively, if a eigenmode oscillates more rapidly in space, more charge transfers across small distances, having lower resistance and hence happening faster. $\hat{\bm{T}}_{a}$ may have very small eigenvalues close to [math] as well, but the magnitudes of the corresponding eigenmodes are also very small. This is because the eigenmodes which rapidly oscillate in space usually don’t correlate with the shape of the total potential distribution. We will show an example of this in Section IV.

III Responses to Typical Stimuli

With the model developed in Section II, we study the system responses to three typical stimuli in electrochemical measurements: sinusoidal waveforms in Section III.1, chronoamperometry (controlled potential method) in Section III.2 and chronopotentiometry (controlled current method) in Section III.3. With the approximation of (21), (14) becomes

[TABLE]

III.1 Sinusoidal Waveforms

A sinusoidal waveform is applied to the electrode: $V(t)=V_{0}e^{j\omega t}$, where $j$ is the imaginary unit. Using the time-domain Fourier transform in (23) yields

[TABLE]

It is not possible to solve (24) explicitly without assuming a specific electrode configuration, but we can estimate the impedance

[TABLE]

at the extremes of the frequency $\omega$.

Now we use $\left\|\cdot\right\|$ to denote the Hilbert–Schmidt norm. As explained in Section II, $\left\|\hat{\bm{T}}_{a}\right\|\gg\left\|\hat{\bm{T}}_{b}\right\|$. In order of magnitude, $\left\|\hat{\bm{R}}_{a}(C_{p}+C_{d})\right\|$ is close to $\left\|\hat{\bm{T}}_{a}\right\|$.

When $\omega\left\|\hat{\bm{T}}_{b}\right\|\gg 1$, $\hat{\bm{Z}}\approx\hat{\bm{R}}_{a}$, and we have

[TABLE]

We see from (26a) that at high frequencies the interface is equipotential, and from (26b) that the current is changing in phase with voltage, and the access resistance associated with the EP boundary condition is measured.

When $\omega\left\|\hat{\bm{T}}_{a}\right\|\ll 1$, $\hat{\bm{Z}}\approx{\left[j\omega(C_{p}+C_{d})\right]}^{-1}$, and we have

[TABLE]

We see from (27a) that at low frequencies, the interface is not equipotential, and from (27b) that the current is shifted by $90^{\circ}$ and is proportional to the total surface capacitance $C_{p}+C_{d}$.

Now if $\omega\left\|\hat{\bm{T}}_{a}\right\|\gg 1$ but $\omega\left\|\hat{\bm{T}}_{b}\right\|\ll 1$, we have

[TABLE]

(28) is more intuitive when $R_{ct}$ and $C_{p}$ are uniform, which gives $\hat{\bm{Z}}=\hat{\bm{R}}_{a}+R_{ct}$. Similar to high frequencies, at middle frequencies, the boundary is equipotential and the current is in phase with voltage. The impedance is the sum of the access resistance and the charge transfer resistance.

III.2 Chronoamperometric Response

To study the transient behavior, we focus on the eigenmodes with long characteristic times, and assume that all eigenmodes of $\hat{\bm{T}}_{b}$ decay infinitely fast. Explicitly, we assume $\hat{\bm{T}}_{b}\partial_{t}+\hat{I}\approx\hat{I}$, and (23) becomes:

[TABLE]

Without loss of generality, we assume $U_{p}(\bm{r},0)=U(\bm{r},0)=0$, and $V(t)=0$ when $t\leq 0$. We take $V(t)$ of the form $V(t)=V_{0}+vt$ for $t>0$. Since $U_{p}(t)$ must be continuous, we have:

[TABLE]

At steady state, we have $\dot{\varphi}=0$. Together, we have:

[TABLE]

Let $\Upsilon_{a,l}$ be the $l^{\mathtt{th}}$ normalized eigenfunction of $\hat{\bm{T}}_{a}$, corresponding to the eigenvalue $\tau_{l}$. We expand $\varphi$ in the basis $\{\Upsilon_{a}\}$, with the coefficients $\{\varphi_{l}\}$. If $C_{p}$ is uniform, $\{\Upsilon_{a}\}$ is orthonormal, and the expansion is straightforward:

[TABLE]

For the more general case when $C_{p}$ is not uniform, we perform the Gram-–Schmidt orthogonalization to $\{\Upsilon_{a}\}$ with a coefficient matrix $\bm{G}$:

[TABLE]

such that $\{\tilde{\Upsilon}_{a}\}$ is orthonormal. We have:

[TABLE]

Note that when $C_{p}$ is uniform, $\bm{G}=\bm{I}$, the identity matrix. By the principle of superposition, we have

[TABLE]

The solution consists of a steady state component $C_{p}(\bm{r})v$, which is PCD, and a transient component consisting of eigenmodes $\{\varphi_{l}\Upsilon_{l}(\bm{r})\}$ that exponentially decay at the rates $\tau_{l}$. We call such transient behavior the EP-PCD transition (or EP-UCD when the interface material is uniform), whose longest characteristic time corresponds to the largest eigenvalue of $\hat{\bm{T}}_{a}$:

[TABLE]

III.3 Chronopotentiometric Response

With controlled total current, (29) still holds, but the solution is not the superposition of exponentially decaying eigenmodes, since $\dot{V}$ is no longer constant in time. Without loss of generality, we assume $U_{p}(\bm{r},0)=U(\bm{r},0)=0$, and the total current $I_{tot}(t)=0$ when $t\leq 0$. We take $I_{tot}(t)=I_{0}$ for $t>0$. Denote $u_{l}$ the net current flow coefficient of the $l^{\mathtt{th}}$ eigenmode:

[TABLE]

so that $\varphi_{l}u_{l}$ is the net current of $\varphi_{l}\Upsilon_{a,l}$. We expand $\hat{\bm{R}}_{a}C_{p}$ in basis $\{\Upsilon_{a}\}$ with coefficients $\{\upsilon_{l}\}$:

[TABLE]

so that $\dot{V}\upsilon_{l}$ is the component of $\hat{\bm{R}}_{a}C_{p}\dot{V}$ in $\Upsilon_{a,l}$.

For the initial and the steady state conditions, (30) and (31) still hold. We combine (31a) and (34), and $V(0^{+})$ is given by $V_{0}$ in

[TABLE]

We can thereby determine $\{\varphi_{l}(0^{+})\}$.

Total current does not change, so we have:

[TABLE]

In basis $\{\Upsilon_{a}\}$, (29) becomes:

[TABLE]

Let $\bm{y}$ be a vector of variables:

[TABLE]

and $\bm{\Gamma}$ be a matrix of coefficients:

[TABLE]

We combine (40) and (41), and have:

[TABLE]

$\Gamma$ is full-rank, and solving for the transient behavior with controlled current becomes a standard problem of homogeneous linear dynamic system, as in (44).

IV Solution for a Disk Electrode

We now consider a disk electrode placed at the center of an insulating plane, with electrolyte filling the half-space above the plane. For simplicity of mathematical forms, from here on, we assume uniform surface material, with constant $R_{ct}$, $C_{p}$ and $C_{d}$. The majority of the theoretical derivation has been done in Nisanciogğlu and Newman (1973a), with the assumption of constant electrochemical reaction potential ($C_{p}\xrightarrow{}\infty$). As a demonstration of our more intuitive framework developed in Sections II and III.2, we apply it to this problem in Section IV.1, and then compare the results with experimental measurements in Section IV.2.

IV.1 Theoretical Derivation

We consider a disk electrode of radius $a$. $R_{P}$, $C_{p}$ and $C_{d}$, as defined in Section II are uniform. In the Cartesian coordinates, we have $E=\{(x,y,z\geq 0)\}$, $A=\{(x,y,0):\sqrt{x^{2}+y^{2}}\leq a\}$ and $D=\{(x,y,0):\sqrt{x^{2}+y^{2}}>a\}$. Following Newman (1966a), we will use the elliptic coordinate system $(\xi,\eta)$. Laplace’s equation in elliptic coordinates is

[TABLE]

with the boundary conditions

[TABLE]

We note that

[TABLE]

By Newman (1966b), the solution to (45) is

[TABLE]

with $l\in\mathbb{N}$, ${k_{l}}$ are constant coefficients, $P_{2l}$ is the $2l^{\mathtt{th}}$ Legendre polynomial of the first kind, and $X_{l}(\xi)$ is the solution to

[TABLE]

We normalize $P_{2l}$, so that $\tilde{P}_{2l}:=\sqrt{4l+1}P_{2l}$ form an orthonormal an orthonormal basis of functions in $\{f\in C^{\infty}:[0,1]\rightarrow\mathbb{R},f^{\prime}(0)=0\}$. $k_{l}$ in (48) are chosen so $\Phi(0,\eta)$ matches a given $\varphi(\eta)$. With $X_{l}(0)=1$, (48) becomes

[TABLE]

By the definition of $\hat{\bm{S}}$ in (11),

[TABLE]

By (47a), we then have

[TABLE]

Equation [18] of Newman (1966b) gave, without derivation, that

[TABLE]

A detailed derivation was provided in Section $6.9$ of Wang (2016). Appendix B provides the derivation in a more rigorous manner, and proves the monotonicity of $X(\xi)$. Intuitively, this shows that potential distribution is monotonic along each hyperbolic line in the elliptic coordinates.

Therefore, the operator $\hat{\bm{S}}$ has the form

[TABLE]

Since $R_{ct}$ and $C_{p}$ are uniform, we now use these two notations as scalars. $\hat{\bm{T}}_{a}$ becomes

[TABLE]

which shares the same eigenspace with $\hat{\bm{S}}$. We define a dimensionless operator

[TABLE]

so that $\hat{\bm{S}}=\dfrac{2}{\pi\rho a}\hat{\bm{A}}$. We have:

[TABLE]

As $|\eta|\leq 1$, the smallest eigenvalue of $\hat{\bm{A}}$ can then be estimated as $\lambda_{\min}\left(\hat{\bm{A}}\right)\geq 1$, and

[TABLE]

The numerically computed first 4 eigenfunctions of $\hat{\bm{A}}$, together with their respective eigenvalues, are shown in FIG. 3. It turns out that the smallest eigenvalue of $\hat{\bm{A}}$ is about $1.8$, which gives

[TABLE]

Note the similarity between $\Phi$ and the electrostatic potential of a charged disk. Specifically, if a flat disk in free space has charge density $\sigma(r)=2\varepsilon_{0}\rho i(r)$, then $\Phi$ is also the potential distribution around the charged disk. With UCD, the potential at the center of the disk electrode is

[TABLE]

Although access resistance is not well defined with UCD since the electrode surface is not equipotential, we can define an effective access resistance as the potential at the center of the electrode divided by the total current:

[TABLE]

This value is higher than the widely accepted access resistance value $\rho/(4a)$ of an equipotential disk electrodeNewman (1966a). There are different definitions of the effective access resistance. For example, Oldham (2004) defines $R_{a,\texttt{eff}}$ as the average of $\varphi$ divided by the total current. We have chosen our definition because the center of the electrode surface is the most typical point to sample when measuring potential in the electrolyte.

To illustrate the dynamics of the current redistribution in chronoamperometry, let $\rho$ and $vC_{p}$, defined in Section III.2, as well as the radius $a$, be unitary, and let $R_{ct}=0$. To keep the total current the same at the initial state and at the steady state, we will choose $V_{0}=\pi\rho\eta Ca/4$. The resulting evolution of the current density and the potential distribution on a disk electrode over time are shown in FIG. 4, and also in the Supplemental Video.

IV.2 Experimental Validation

To experimentally verify whether the dynamics of the total current on a capacitive electrode-electrolyte interface matches the solution described by (35), we performed chronoamperometric measurements. From here on, all potentials are referred to the Ag/AgCl electrode, unless noted otherwise. In order to sustain higher current within the relatively low voltage window, we used an electrode coated with a sputtered iridium oxide film (SIROF) – a material known for its large charge injection capacity (CIC)Cogan et al. (2009). The continuous iridium valency of SIROF between [math] to $0.8\text{\,}\mathrm{V}\text{/}$Pauporté et al. (1999), together with its porous surfaceCogan (2008), enables a large capacitance.

We used an electrochemical cell of the 3-electrode configuration. The working electrode is a $80\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$-diameter platinum disk coated with $400\text{\,}\mathrm{nm}\text{/}$ of SIROF. The electrode was treated with $4$% NaClO solution and plasma cleaning, following the protocol of Duru (2019) (Sections $2.15$ and $3.3$). A large ($>$1\text{,}{\mathrm{cm}}^{2}\text{/}) platinum grid was used as the counter electrode. The reference was an Ag/AgCl electrode in $3$\mathrm{M} KCl solution. The electrolyte is $6$-time diluted phosphate buffered saline (PBS) solution, whose resistivity is $353\text{\,}\mathrm{\SIUnitSymbolOhm}\cdot\mathrm{c}\mathrm{m}$, measured with an electrical conductivity meter.

First, we validated that the electrode kinetics is invariant within the potential range, and that concentration polarization is not the dominating factors in the electrode impedance. As shown in FIG. 5, the black solid line in the top panel represents a step voltage pulse. The corresponding total current is shown by the black solid line of the bottom panel (Trial 1). To show that the electrode kinetics is not varying with potential, we offset the voltage pulse up and down by $100\text{\,}\mathrm{mV}\text{/}$ (Trials 2 and 3), and observed that the current did not change. To check whether the concentration polarization affects the current amplitude, we scaled the voltage pulse by a factor of 2. The current nearly doubled as well, indicating that concentration polarization is negligible since the concentration overpotential does not scale linearly with the current density.

We observed very small charge transfer resistance on this electrode, which allows to neglect the effect of $R_{ct}$ and combine $C_{p}$ and $C_{d}$ into one supercapacitance $C_{s}$. To confirm, we performed electrochemical impedance spectroscopy (EIS), and used the Levenberg-Marquardt method to fit the Bode plot to the circuit diagram in FIG. 6(a). The measurement and the fitting curves are plotted in FIG. 6(b). From the fitting, $R_{ct}+R_{a}=$22.3\text{,}\mathrm{k\SIUnitSymbolOhm}\text{/}$, within$2$\%$ error range of the EP access resistance predicted by $\rho/(4a)=$22.0\text{,}\mathrm{k\SIUnitSymbolOhm}. Per the discussion in Section [III.1](#S3.SS1), we confirm that $R_{ct}$ is negligible. We also found $C_{s}=$8.52\text{\,}\mathrm{mF}\text{/}{\mathrm{cm}}^{2} from the fitting.

A voltage waveform including a step and a ramp, with $V_{0}=$-100\text{,}\mathrm{mV}\text{/} and $v=$-3.24\text{\,}\mathrm{V}\text{/}\mathrm{s}, as defined in Section III.2, was applied to the SIROF electrode (top panel in FIG. 4). Note that this waveform is different from the one used in FIG. 4, and hence the total current at the initial and the steady states is not the same either. The resulting current waveform is shown in the same panel.

One competing theory is that the surface instantly exchanges charge laterally, so the charge accumulation in the capacitor is uniform and the interface is always EP. This is the assumption behind the RC fitting in Myland and Oldham (2005); Behrend et al. (2008); Boinagrov et al. (2015). To compare the experimental results with predictions of the constant EP theory and of the EP-UCD transition, time derivative of the total current was calculated and plotted for both models and the measurement. The EP-UCD perdiction is from (35), while the constant EP perdiction is a simple RC process. As can be seen in the bottom panel of FIG. 7, the measurement matches our theory rather than the constant EP assumption.

V Discussion

A highly conductive electrode is always equipotential in its bulk, but this equipotentiality is often confused with the surface layer of electrolyte at the interface, which is the boundary typically modeled as the Helmholtz plane. Electrode kinetics and concentration polarization are the two mechanisms previously considered to cause uniform secondary current distribution. Another important mechanism is the charge accumulation on the interface, which is often under-appreciated in applications since it is not reflected in the initial current distribution. Previously, this effect has only been modeled for disk electrodes of uniform surfaces, assuming constant reaction potentials (Nisanciogğlu and Newman (1973a, b)) or no electrochemical reactions (Myland and Oldham (2005); Behrend et al. (2008)).

As we show in Section II, for any geometry and any combination of surface materials, the current density eventually reaches the PCD steady state (or UCD if the capacitance per unit area is the same over the whole electrode). Redistribution of the current from the initial non-uniform spread at an equipotential state is driven by the uneven charge accumulation at the capacitive interface until it reaches the PCD, when the potential of all parts of the interface rises at the same rate. This transition is described by the superposition of exponentially decaying eigenmodes, each of which has a different time constant. Each eigenmode is a surface potential distribution that elicits the circuit response to change itself proportionally. The shorter the time constant is, the faster the eigenmode decays. For a disk electrode of radius $a$ and with uniform surface capacitance $C$, the dominant (longest) time constant is $0.864\rho Ca$, only $10\%$ larger than the simple RC time constant $\tau_{EP}=\pi\rho Ca/4=0.785\rho Ca$, where the EP access resistance $R_{EP}=\rho/(4a)$ is assumed. In an earlier finite element modelingMyland and Oldham (2005), the total current was fit to one simple RC process, resulting in $8.7\%$ increment of the time constant compared to $\tau_{EP}$, which roughly matches our result. Strictly speaking, there are two sets of eigenmodes and time constants, dominated by the pseudocapacitance $C_{p}$ and the double-layer capacitance $C_{d}$, respectively. However, since the latter is faster than the former, when studying the transient behavior, we consider the latter to be instant with negligible effect on the circuit behavior. This requires $C_{p}\gg C_{d}$ everywhere, which may not be true if a surface consists of both electrochemically active parts and inert parts, but we can avoid this subtlety by choosing, nominally, $C_{d}=0$, $R_{ct}=0$ and $C_{p}$ the double-layer capacitance. Empirically, the inverse time constants of different eigenmodes are separated almost evenly, which is equivalent to an asymptotic approximation conjectured by Troesch and TroeschTroesch and Troesch (1972) in a solution of a problem in fluid dynamics, and confirmed computationally up to the $200^{\mathtt{th}}$ time constantWang (2016).

It is important to note that the measurements of the access resistance using electrical impedance spectroscopy (EIS) correspond to the high end of the frequency range. At frequencies exceeding the inverse time constants of the current redistribution, the interface remains practically equipotential. Therefore, the access resistance measured in EIS is associated with the EP boundary condition. In the middle of the frequency range, we should see the sum of the EP acess resistnace and the charge transfer resistance, if $C_{p}$ and $C_{d}$ are separated sufficiently apart.

Under the constant EP boundary condition, the current in response to a voltage step with a ramp is a simple exponential decay to the steady state, with a time constant of RC. Distinguishing this curve from a plot corresponding to our theory is not easy since the dominant time constant is only slightly longer than in the constant EP theory, and the magnitude of the slowest eigenmode is the largest. Therefore, the bottom panel in Figure 7 compares the time derivative of the total current. Since the constant EP theory has only one decaying mode, while the EP-UCD transition has infinitely many and much faster decaying eigenmodes, total current decreases faster at the beginning of the pulse, as can be seen in the plot.

At steady state, the PCD boundary condition enables control of the current distribution on various parts of the interface by selecting electrode materials of different capacitance per unit area. For example, if a part of the Au electrode is coated with SIROF, and the pulse duration exceeds the characteristic EP-PCD transition time, the current will flow primarily through the SIROF area, while the Au surface will be practically passive since its capacitance is about $1000$ times smaller than that of SIROF. This effect was discovered in Flores et al. (2016) but only analyzed using a discrete circuit approximation. The phenomenon of PCD greatly simplifies the 3-D electrode fabrication by electroplating: the side walls of the Au-electroplated electrode do not have to be coated with an insulator. They can remain exposed to the liquid since the SIROF on top of these walls will collect vast majority of the currentFlores et al. (2018a, b); Ho et al. (2018). Similarly, leads to a high-capacitance electrode do not have to be well-insulated from the medium as long as their capacitance is much smaller than that of the target electrode. For example, the electrodes used in Musk et al. (2019).

Understanding the distribution of electric field in the medium is particularly important for proper design of the electro-neural interfaces. For example, if the pulse duration is significantly shorter than the EP-PCD transition time, the electric current will flow primarily from the electrode edges. This will result in highly enhanced electric field in these areas, which may stimulate and even damage the nearby cells much more than the average current density calculated by dividing the total current by the total electrode areaWang et al. (2014). The edge effect can be effectively avoided if the electrode capacitance is selected such that the characteristic transition time is below the intended pulse duration. In addition, the electrode capacitance can be gradually reduced toward the edges, for example, by decreasing the SIROF thickness using partial shadowing techniques.

VI Conclusions

We provided an analytical solution describing the dynamics of the current redistribution on capacitive electrode-electrolyte interfaces and validated our theory experimentally. We demonstrated that current and voltage redistribute over time from the initial non-uniform spread to the steady state, where the current density at the surface is proportional to the capacitance per unit area. This transition can be described as a superposition of the exponentially decaying eigenmodes. The slowest and dominant eigenmode of a disk electrode has a time constant similar to RC of the electrode. We also note that since the EIS based measurements of the access resistance are performed at high frequencies, they correspond to equipotential boundary condition, which is different from the access resistance at low frequencies. To avoid the strong edge effects on large electrodes, the capacitance of the electrode material should be selected so that the EP-PCD transition time does not significantly exceed the intended pulse duration.

Acknowledgements.

Authors would like to thank Prof. Christopher Chidsey from Stanford University and Dr. Boshuo Wang from Duke University for the very helpful discussions. Authors would also like to thank Yibin Shu from the Institute of Education at Tsinghua University in China for his help with preparation of the key image. Funding was provided by the National Institutes of Health (Grants R01-EY-018608, R01-EY-027786), Stanford Neurosciences Institute, and Research to Prevent Blindness.

Appendix A

Here we show that the operator $\hat{\bm{S}}$ defined in (11) is positive-definite.

Proof.

Let $\Phi$ and $\Psi$ be two non-zero potential distributions in $E$ and define

[TABLE]

with similar definitions for $\psi_{0}$ and $\psi_{1}$. By (4) and (5)

[TABLE]

We may now write

[TABLE]

The last equality above used (63a). By the divergence theorem and (10), we now have

[TABLE]

thus

[TABLE]

It follows that $\hat{\bm{S}}$ is Hermitian, which is the result of the Lorentz reciprocity. Furthermore, as

[TABLE]

we have

[TABLE]

and it follows that $\hat{\bm{S}}$ is positive-definite.

∎

Appendix B

The general solution to (49) is

[TABLE]

where $j=\sqrt{-1}$, $c_{1},c_{2}\in\mathbb{C}$ are coefficients and $Q_{2l}$ is the $2l^{\mathtt{th}}$ Legendre polynomial of the second kind. The boundary conditions are

[TABLE]

As $Q_{2l}(0)=0$, by (70a) we have $c_{1}=\frac{1}{P_{2l}(0)}$. By equation (12.216) in Arfken and Weber (2005), we also have

[TABLE]

where $R_{2l-1}$ is a polynomial of degree $(2l-1)$ with only odd-order terms, thus

[TABLE]

Using (70b), we conclude that

[TABLE]

and

[TABLE]

By Theorem 1 below, a solution satisfying (70) exists and it must be of the form

[TABLE]

Now because $P^{\prime}_{2l}(0)=0$, it follows that

[TABLE]

Theorem 1**.**

A monotonically decreasing solution to (49) satisfying the boundary conditions (70) exists.

Proof.

When $l=0$, $X_{l}(\xi)=1-\frac{2}{\pi}\arctan\xi$ is a valid solution. so we assume $l\geq 1$.

Let $X_{n}$ be a sequence of solutions to (49), defined on $0\leq\xi\leq n$, with the boundary conditions

[TABLE]

Choosing $c_{1}={1}/{P_{2l}(0)}$ and $c_{2}=-{P_{2l}(in)}/{Q_{2l}(in)}$ gives an explicit form for $X_{n}$.

We claim that $X_{n}$ monotonically decreases in $(0,n)$, and prove this by contradiction.

Indeed, (49) gives

[TABLE]

If $X_{n}$ is not monotonic, there exists a local extremum $\xi_{0}\in(0,n)$ such that $X_{n}^{\prime}(\xi_{0})=0$. We see from (78) that $X_{n}^{\prime\prime}(\xi_{0})$ has the same sign as $X_{n}(\xi)$. Thus, if $X_{n}(\xi_{0})>0$, then $\xi_{0}$ is a local minimum. Therefore, there exists a local maximum $\xi_{1}\in(\xi_{0},n)$ such that $X_{n}(\xi_{1})>0$. However, at $\xi_{1}$, which is also an extremum, $X_{n}^{\prime\prime}(\xi_{1})$ and $X_{n}(\xi)$ have different signs, which is a contradiction. Similarly, assuming $X_{n}(\xi_{0})<0$ leads to a contradiction as well. Therefore, $X_{n}$ monotonically decreases in $(0,n)$.

We also claim that $X_{n}(\xi)$ is increasing in $n$, that is, if $m>n$ then $X_{m}(\xi)\geq X_{n}(\xi)$ for $0\leq\xi\leq n$. Indeed, monotonicity of $X_{m}(\xi)$ in $\xi$ implies that $X_{m}(n)>0$, thus $Z=X_{m}(\xi)-X_{n}(\xi)$ satisfies (49) with $Z(0)=0$ and $Z(n)>0$. By the same argument, $Z$ can not attain a negative minimum, thus $Z(\xi)>0$ for all $0<\xi<n$. Therefore, the limit

[TABLE]

exists, and $X_{l}\geq 0$ is monotonically decreasing. We claim that

[TABLE]

and prove this by contradiction.

Assume there exists $\epsilon>0$ such that $X_{l}(\xi)\geq\epsilon$ for all $\xi\geq 0$. The integral of the left side of (49) yields

[TABLE]

The integral of the right side of (49) yields

[TABLE]

As $t\rightarrow+\infty$, we have $2l(l+1)\epsilon t>-X^{\prime}_{l}(0)$, which is a contradiction, and (80) follows. ∎