Also at ]Center for Advanced Materials Research, Research Institute of
Sciences and Engineering, University of Sharjah, Sharjah,, P.O. Box 27272,
United Arab Emirates

# Modified Poisson-Nernst-Planck theory for low-to-mid frequency immittance of
electric double-layer capacitors

Anis Allagui∗<EMAIL_ADDRESS>Dept. of Sustainable and Renewable Energy
Engineering, University of Sharjah, Sharjah, P.O. Box 27272, United Arab
Emirates  [  Dept. of Mechanical and Materials Engineering, Florida
International University, Miami, FL33174, United States  Hachemi Benaoum
Dept. of Applied Physics and Astronomy, University of Sharjah, PO Box 27272,
Sharjah, United Arab Emirates  Hichem Eleuch  Dept. of Applied Physics and
Astronomy, University of Sharjah, PO Box 27272, Sharjah, United Arab Emirates
Chunlei Wang  Dept. of Mechanical and Materials Engineering, Florida
International University, Miami, FL33174, United States

###### Abstract

Understanding the system-level spectral immittance response of capacitive
energy storage devices with analytically tractable physics-based models is not
only important for the progress of the technology, but also allows to develop
new physical insights more easily. Here, we report a modified
Poisson–Nernst–Planck (PNP) system describing charge concentration and
electric potential as a model of electrokinetics for electrodes showing mixed
resistive-capacitive behavior. This is done by (i) incorporating time shifts
between the current fluxes and both concentration gradients of charged species
and the electric field, and (ii) introducing time fractional derivatives in
the continuity equation. The aim is to characterize the deviation of
immittance from that of ideal capacitors both at close-to-dc frequencies where
the impedance angle for example is larger than -90 deg., and also at mid-range
frequencies where the system veers progressively toward resistive behavior.
This latter tendency is important to model in order to identify the extend of
the capacitive bandwidth of the device from the rest. Solution and simulation
results to the one-dimensional modified PNP system for symmetric
electrolyte/blocking electrode configuration are presented and discussed.

##  I Introduction

The analysis and interpretation of measured spectral immittance profile is
commonly carried out using equivalent circuit models. In particular, the non-
ideal behavior that cannot be properly captured by circuits of $R$, $L$ and
$C$ elements is usually described by circuits containing the constant phase
element (CPE) Allagui  _et al._ (2020a, 2016a, 2018a, 2016b); Allagui,
Elwakil, and Psychalinos (2021). The CPE, also known as mono-order fractional
capacitor, is widely employed for the modeling of electric double-layer
capacitors (EDLC) and porous electrodes Bisquert  _et al._ (1998); Allagui
_et al._ (2016a, 2018a); Allagui, Elwakil, and Fouda (2021); Allagui and Fouda
(2021); Allagui  _et al._ (2020b); Allagui, Zhang, and Elwakil (2018);
Baboukani  _et al._ (2019); Khakpour  _et al._ (2019), batteries Allagui  _et
al._ (2021); Cheng  _et al._ (2018); Xiong  _et al._ (2018); Cuervo-Reyes  _et
al._ (2015), solar cells Zhang  _et al._ (2019, 2020), and sensor devices
Borini  _et al._ (2013); Song  _et al._ (2020); Meunier  _et al._ (2020). It
has the (empirical) power-law complex impedance function of the form:

$Z(s)=V(s)/I(s)=1/(s^{\alpha}C_{\alpha})$ (1)

where $C_{\alpha}$ is in units of F sα-1 (pseudocapacitance), $0<\alpha<1$ and
$s=j\omega$ the complex angular frequency. The real and imaginary parts of a
CPE are $\mathrm{Re}(Z)=\omega^{-\alpha}C_{\alpha}^{-1}\cos(\alpha\pi/2)$, and
$\mathrm{Im}(Z)=-\omega^{-\alpha}C_{\alpha}^{-1}\sin(\alpha\pi/2)$,
respectively, and thus a constant, frequency-independent impedance phase angle
$\phi(Z)=-\alpha\pi/2$, which can take values anywhere between -90 deg. and
zero. We note that in the time-domain, the CPE is characterized by the
current-voltage fractional order differential equation Allagui  _et al._
(2018a, 2016a, 2020b); Allagui and Elwakil (2021):

$i(t)=C_{\alpha}\,{}_{0}D_{t}^{\alpha}v(t)$ (2)

where ${}_{0}D_{t}^{\alpha}v(t)$ is defined here in the Caputo sense by
Podlubny (1998):

${}_{0}D_{t}^{\alpha}v(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{\mathrm{d}v(\tau)/\mathrm{d}\tau}{(t-\tau)^{\alpha}}d\tau$
(3)

For the limiting case of $\alpha=0$, the impedance in Eq. 1 is real-valued,
$Z=R$ (i.e. $1/C_{\alpha}\leftarrow R$) for any frequency, and Eq. 3 turns to
be Ohm’s law $i(t)=v(t)/R$ (zero-order derivative ${}_{0}{D}_{t}^{0}f(t)$
gives $f(t)$, i.e. the identity operator). Whereas for $\alpha=1$, the
impedance is that of an ideal capacitor $1/(sC_{1})$ and the current-voltage
characteristic equation is $i(t)=C_{1}\mathrm{d}v(t)/\mathrm{d}t$. The Warburg
impedance is in fact a special case of the CPE corresponding to $\alpha=1/2$.
If the value of $\alpha$ is taken from $\left]-1;0\right[$, we may speak of
fractional inductors of mixed behavior lying between that of ideal inductors
and resistors Fouda  _et al._ (2016), but that is not the focus of this work.

Figure 1: Results of open-circuit electrochemical impedance spectroscopy of an
all-solid-state (symmetric) multi-layer graphene electric double-layer
capacitor Allagui  _et al._ (2018b, 2017): (a) Nyquist plot of imaginary part
vs. real part (frequency is swept as a parameter), and (b) impedance phase
angle vs. log of frequency.

Because of the extra degree of freedom in the CPE impedance, that is the
dispersion coefficient $\alpha$, this element is proved successful for
describing phenomenologically the inclined and depressed Nyquist impedance
plots of real vs. imaginary parts of many systems, as well as their time-
domain responses Allagui  _et al._ (2016b, 2021); Allagui, Elwakil, and
Psychalinos (2021); Allagui  _et al._ (2016a); Fouda  _et al._ (2016); Zhang
_et al._ (2019). In fact the pure capacitive behavior is rarely encountered in
real measurements which made the CPE a very successful fitting tool. However,
typical experimental results, as shown in Fig. 1 for multi-layered graphene-
based EDLC, exhibit fractional order behavior in the close-to-dc region
(impedance phase angle is practically constant and equal to -83∘ over the
frequency bandwidth 0.01 to $<1$ Hz), but as the frequency is increased we
observe an increase of resistive behavior at the expense of the capacitive
behavior. The low-frequency segment can be fitted well enough with a single
CPE or a CPE with a resistor, but that cannot be the case anymore as the
frequency is increased. It is usually required to have additional elements in
the modeling circuit which can complicate the analysis and hinder the physical
tractability to the system under consideraion. Furthermore, the CPE electrical
parameters have still unsatisfactory interpretations, and their relationship
to real physical measurables is still a subject of debate. The current
understanding revolves around subdiffusive charge transport where the mean
square displacement of a particle follows a power law function of time
Kosztołowicz and Lewandowska (2009); Allagui  _et al._ (2020a), fractal
morphologies and interfacial roughness of electrodes Kant, Kumar, and Yadav
(2008); Maritan, Stella, and Toigo (1989), or some form of distribution of
response relaxation Zhang  _et al._ (2015).

The motivation behind this work is to study the anomalous frequency-domain
immittance function observed with EDLCs at a system level using physics-based
models. Our goal is to describe both the deviation of the frequency response
of such devices from that of ideal capacitors at low frequencies (i.e. CPE
behavior), but also to show the evolution of the system towards resistive
behavior as the frequency is increased. To this end we employ the standard
electrokinetic model based on Poisson-Nernst-Planck (PNP) equations, which is
well-known for treating electrochemical transport of charged species by
relating their fluxes to both concentration gradients and the electric field
Barbero (2017). It has been applied and validated in electrochemistry and
analytical chemistry, sensor devices, and in biology Kilic, Bazant, and Ajdari
(2007); Schmuck and Bazant (2015); Alexe-Ionescu  _et al._ (2020); Koklu,
Ajaev, and Beskok (2019); Barbero (2017). In addition, these equations also
appear in plasma physics, and semiconductor device modeling, where they are
known as the diffusion-drift equations or the van Roosbroeck equations. We
consider the case of perfect blocking electrodes without convective
electrolyte transport, steric effects due to finite ion sizes, or Faradaic
electrochemical reactions at the electrodes. We will, however, incorporate a
time shift between the fluxes of charged species and both their concentration
gradients and the electric field in a similar way as Cattaneo’s constitutive
formulation Kosztołowicz and Lewandowska (2009); Compte and Metzler (1997). We
recall that in Cattaneo’s formulation of Fickian diffusion for instance, the
constitutive equation relating flux $\mathbf{J}(\mathbf{x},t)$ to the
distribution function of the diffusing species $c(\mathbf{x},t)$, i.e.

$\mathbf{J}(\mathbf{x},t)=-D\,\nabla c(\mathbf{x},t)$ (4)

would be considered (up to the first order) as:

$\mathbf{J}(\mathbf{x},t)+\tau_{c}\frac{\partial\mathbf{J}(\mathbf{x},t)}{\partial
t}=-D\,\nabla\mathbf{c}(\mathbf{x},t)$ (5)

with non-zero value of the relaxation time constant $\tau_{c}$ (typically in
the range of $10\sim 100$ s for porous materials and biological tissues,
$10^{-10}\sim 10^{-14}$ s for metals, $10^{-8}\sim 10^{-10}$ s for gases). The
phenomenological derivation of this equation is based on the assumption that
the flux of particles is not generated by the concentration gradient
instantaneously at time $t$, as it is the case with the normal diffusion
equation, but is delayed by a certain time lag $\tau_{c}$ Lewandowska and
Kosztołowicz (2008). Such an approach has been extensively studied and
verified in heat transfer (Fourier’s law) Ruggeri, Muracchini, and Seccia
(1990); Joseph and Preziosi (1989) and diffusion problems (Fick’s law) Compte
and Metzler (1997); Kosztołowicz (2014); Kosztołowicz and Lewandowska (2009);
Garcia-Colin and Olivares-Robles (1995); Qi and Jiang (2011). Without this
modification, the fundamental solution to the (parabolic) diffusion equation
(obtained by substituting Eq. 4 in the continuity equation, ${\partial
c_{i}(\mathbf{x},t)}/{\partial t}+\nabla\mathbf{J}_{i}(\mathbf{x},t)=0$, Eq.
12 below), i.e.

$\frac{\partial c(\mathbf{x},t)}{\partial t}=D\nabla^{2}c(\mathbf{x},t)$ (6)

in response to an initial Dirac delta function
$c(\mathbf{x},t=0)=\delta(\mathbf{x})$, is given by the Gaussian distribution:

$c(\mathbf{x},t)=(4\pi Dt)^{-d/2}e^{-\textbf{x}^{2}/(4Dt)}$ (7)

Eq. 7 clearly leads to non-zero, finite amount of diffusing species at large
distance away from the origin even for very small values of $t$ Compte and
Metzler (1997). This can be interpreted as a propagation of species taking
place at infinite velocity, which is unphysical. Note that combining Eq. 5
with the unchanged continuity equation ${\partial c(\mathbf{x},t)}/{\partial
t}=-D\nabla\mathbf{J}(\mathbf{x},t)$, one obtains the second-order hyperbolic
differential equation of diffusion as per Cattaneo’s formulation Compte and
Metzler (1997); Kosztołowicz and Lewandowska (2009):

$\frac{\partial c(\mathbf{x},t)}{\partial
t}+\tau_{c}\frac{\partial^{2}c(\mathbf{x},t)}{\partial
t^{2}}=D\nabla^{2}c(\mathbf{x},t)$ (8)

In the PNP system we propose below for describing the frequency behavior of
EDLCs, we will adopt a similar approach which will provide, as we shall show,
the tools for capturing the tendency of the electrochemical capacitor towards
resistive behavior as the frequency is increased. Furthermore, we consider
fractional-order time derivative in the continuity, which is commonly used for
describing the anomalous behavior of complex system when integer-order models
fails, such as in disordered systems, gels and porous electrodes Kosztołowicz
and Lewandowska (2009); Gomez-Zamudio  _et al._ (2016); Lenzi  _et al._
(2021); Qi and Jiang (2011). Such a generalization is connected to Hilfer’s
fractional invariance or fractional stationarity Hilfer (1995).

The remainder of the article is organized as follows. In Section II, we
describe the classical ion transport in an electrochemical system based on PNP
theory, and the associated boundary conditions we considered to obtain an
expression for the admittance/impedance functions of ideal capacitors. We also
derive the new expressions for the admittance/impedance functions (Eq. 36) by
solving the modified PNP system. Simulations and discussion of the results are
presented in Section III. Finally, in Section IV, we present our concluding
remarks highlighting the main advantages of the modified PNP model for
describing anomalous impedance of fractional EDLCs over extended frequency
bandwidths.

##  II Theory

###  II.1 Classical and modified PNP systems

First, we present the basic equations used for describing electrified
electrodes in contact with electrolytes, which are the flux equations, the
continuity equation and the local electroneutrality assumption Manzanares and
Kontturi (2003). Assuming constant pressure and temperature, the flux equation
for a given species $i$ (of positive or negative charge) in the electrolyte is
derived from the difference in electrochemical potentials
$\mu^{\prime}_{i}-\mu^{\prime\prime}_{i}$ between two spatially adjacent
regions of volume elements that we denote by $V^{\prime}$ and
$V^{\prime\prime}$. This is related to the change in internal energy of these
two volumes which are
$\text{d}U^{\prime}=T\text{d}S^{\prime}+\sum_{i}\mu_{i}\text{d}n^{\prime}_{i}$
and
$\text{d}U^{\prime\prime}=T\text{d}S^{\prime\prime}+\sum_{i}\mu_{i}\text{d}n^{\prime\prime}_{i}$,
with the additional assumption of transport without energy exchange with the
surroundings, i.e., $\text{d}U^{\prime}+\text{d}U^{\prime\prime}=0$ and
$\text{d}n^{\prime}_{i}+\text{d}n^{\prime\prime}_{i}=0$. Thus, from the Second
Law of Thermodynamics for the control volume $V^{\prime}\cup
V^{\prime\prime}$, we have
$T(\text{d}S^{\prime}+\text{d}S^{\prime\prime})=\sum_{i}(\mu^{\prime\prime}_{i}-\mu^{\prime}_{i})\text{d}n^{\prime}_{i}\geqslant
0$, which indicates that $\text{d}n^{\prime}_{i}$ is determined only by the
difference $\mu^{\prime\prime}_{i}-\mu^{\prime}_{i}$ when the transport of
species $i$ is not coupled to that of other species. If we denote by
$\mathbf{n}\,\text{d}A$ the area separating the two volumes and oriented along
the unit vector $\mathbf{n}$ from $V^{\prime}$ and $V^{\prime\prime}$, then
the number of moles of species $i$ crossing the surface in a time $\text{d}t$
is $\text{d}n^{\prime}_{i}=\
\mathbf{J}_{i}\,\mathbf{n}\,\text{d}A\,\text{d}t$. When the difference in
electrochemical potentials is not too large, we can assume that
$\text{d}n^{\prime}_{i}/\text{d}t$ is proportional to the gradient of the
electrochemical potential normal to the surface (i.e.
$\nabla\mu_{i}\,\mathbf{n}$), and thus express the vector flux density as:

$\mathbf{J}_{i}(\mathbf{x},t)=c_{i}(\mathbf{x},t)\mathbf{v}_{i}(\mathbf{x},t)=-\frac{D_{i}c_{i}(\mathbf{x},t)}{RT}\nabla\mu_{i}(\mathbf{x},t)$
(9)

Here $\mathbf{v}_{i}(\mathbf{x},t)=-u_{i}\nabla\mu_{i}(\mathbf{x},t)$ is the
velocity of species $i$ of composition $c_{i}$, where $u_{i}=D_{i}/RT$
(Einstein relation) is its mobility and $D_{i}$ its diffusion coefficient
(thermal motion). Furthermore, $\nabla\mu_{i}(\mathbf{x},t)$ can be related to
the changes in composition $c_{i}(\mathbf{x},t)$ and in electric potential
$\phi(\mathbf{x},t)$ such that:

$\nabla\mu_{i}(\mathbf{x},t)=RT\nabla\ln
c_{i}(\mathbf{x},t)+z_{i}F\nabla\phi(\mathbf{x},t)$ (10)

where $F$ is the Faraday constant and $z_{i}$ is the charge number of species
$i$. This leads to the standard Nernst-Plank ionic flux equation of diffusion
and migration in the domain (without convective contribution or Faradic charge
transfer) being:

$\mathbf{J}_{i}(\mathbf{x},t)=-D_{i}\left(\nabla
c_{i}(\mathbf{x},t)+\frac{z_{i}F}{RT}c_{i}\nabla\phi(\mathbf{x},t)\right)$
(11)

To complete the necessary system of equations for describing the
electrochemical system, the time-dependence of composition is given by the
divergence of the flux (continuity equation) such that:

$\frac{\partial c_{i}(\mathbf{x},t)}{\partial
t}+\nabla\mathbf{J}_{i}(\mathbf{x},t)=0$ (12)

with the assumption that no chemical reactions can take place. In addition,
the electrostatic potential distribution in the domain can be determined by
Poisson’s equation as:

$\nabla^{2}\phi(\mathbf{x},t)=-\frac{\rho(\mathbf{x},t)}{\epsilon}=-\frac{F}{\epsilon}\sum_{i}z_{i}c_{i}(\mathbf{x},t)$
(13)

where $\rho(\mathbf{x},t)=\sum_{i}z_{i}Fc_{i}(\mathbf{x},t)$ is the volume
charge density, and $\epsilon$ is the dielectric permittivity of the bulk
solution which is assumed to be constant and frequency-independent. Here the
electric field is considered conservative and related to the potential via
$\mathbf{E}(\mathbf{x},t)=-\nabla\phi(\mathbf{x},t)$.

At this point we introduce the modified PNP system that we propose taking into
account Cattaneo’s approach. In order to consider a finite velocity of
propagation of diffusing species and relax the flux, Cattaneo proposed to
replace the standard linear response (Eq. 4) by Eq. 5 to obtain Eq. 8.
However, the left-hand side of Eq. 5 is nothing but the first-order
approximation of the delayed flux $\mathbf{J}(\mathbf{x},t+\tau_{c})$ in the
limit of small values of $\tau_{c}$, i.e. Kosztołowicz and Lewandowska (2009);
Lewandowska and Kosztołowicz (2008):

$\displaystyle\mathbf{J}(\mathbf{x},t+\tau_{c})$ $\displaystyle=$
$\displaystyle e^{\tau_{c}\frac{\partial}{\partial
t}}\mathbb{J}_{i}\left(\mathbb{x},t\right)$ (14) $\displaystyle\approx$
$\displaystyle\mathbf{J}(\mathbf{x},t)+\tau_{c}\frac{\partial\mathbf{J}(\mathbf{x},t)}{\partial
t}+\ldots$

from which we write:

$\mathbf{J}_{i_{\alpha}}(\mathbf{x},t+\tau_{c})=-\frac{{D_{i_{\alpha}}}c_{i}(\mathbf{x},t)}{RT}\nabla\mu_{i}(\mathbf{x},t)$
(15)

or

$\mathbf{J}_{i_{\alpha}}(\mathbf{x},t)=-\frac{{D_{i_{\alpha}}}c_{i}(\mathbf{x},t-\tau_{c})}{RT}\nabla\mu_{i}(\mathbf{x},t-\tau_{c})$
(16)

Combined with the fractional continuity equation Compte and Metzler (1997);
Hilfer (1995); Qi and Jiang (2011); Hilfer (2000); Kosztołowicz and
Lewandowska (2009),

${}_{0}D_{t}^{\alpha}c_{i}{(\mathbf{x},t)}=-\nabla\mathbf{J}_{i_{\alpha}}(\mathbf{x},t)\quad(0<\alpha<1)$
(17)

the proposed modified PNP system is written in terms of a single lag
$\tau_{c}$ and fractional-order time derivative as follows:

${}_{0}D_{t}^{\alpha}c_{i}{(\mathbf{x},t)}=-\nabla\left\\{\frac{{D_{i_{\alpha}}}c_{i}(\mathbf{x},t-\tau_{c})}{RT}\nabla\left\\{RT\ln
c_{i}(\mathbf{x},t-\tau_{c})+{z_{i}F}\phi(\mathbf{x},t-\tau_{c})\right\\}\right\\}$
(18)

or in terms of two phase lags $\tau_{1}$ and $\tau_{2}$ (for mathematical
convenience, see Eq. 29 below) as Tzou (1995); Compte and Metzler (1997):

${}_{0}D_{t}^{\alpha}c_{i}{(\mathbf{x},t+\tau_{1})}=-\nabla\left\\{\frac{{D_{i_{\alpha}}}c_{i}(\mathbf{x},t+\tau_{2})}{RT}\nabla\left\\{RT\ln
c_{i}(\mathbf{x},t+\tau_{2})+{z_{i}F}\phi(\mathbf{x},t+\tau_{2})\right\\}\right\\}$
(19)

together with the (unchanged) Poisson equation:

$\nabla^{2}\phi(\mathbf{x},t)=-\frac{\rho(\mathbf{x},t)}{\epsilon}$ (20)

The generalization of the continuity equation, that now includes fractional
derivative of constant order $\alpha\neq 1$, stands in direct connection with
Hilfer’s fractional invariance or fractional stationarity Hilfer (1995). It
has been also adopted in various forms, for instance by Compte and Metzler
Compte and Metzler (1997) or Górska Górska (2021) for the study of anomalous
diffusive transport within the framework of Cattaneo, and similarly by
Kosztołowicz and Lewandowska Kosztołowicz and Lewandowska (2009) and others
Lenzi, Evangelista, and Barbero (2009) for study of subdiffusive impedance
spectra, by Kosztołowicz Kosztołowicz (2014, 2019) for the study of Cattaneo-
type subdiffusion-reaction and anomalous diffusion-absorption processes. Note
that in order to maintain proper dimensions of both sides of Eq. 19, the
modified diffusion coefficients $D_{i_{\alpha}}$ are in units of m2/sα.
Furthermore, it is assumed with the suggested model that there is delayed
flux-force relation between the current flux on the one hand and changes in
concentration of species $i$ and in the electric potential in the other hand,
which is taken into account via the net time shift
$\Delta\tau=\tau_{c}=\tau_{2}-\tau_{1}$. This can be particularly justified
for the case of spectral impedance of electrochemical systems, especially at
high frequencies at which the concentrations of species are expected to
oscillate very rapidly Lewandowska and Kosztołowicz (2008). Excitation and
response cannot happen exactly at the same time. More details on the time
shift $\Delta\tau$ will be discussed in Section. III. Furthermore, it is clear
from Catteneo’s formalism (Eq. 5) that the delayed concentration is considered
up to the first order only, whereas in our generalized formulation we take
into account all orders without approximation. Eq. 19 simplifies to
Cattaneo’s-like formulation when taking first order only in terms of
$\tau_{1}$ and setting $\tau_{2}=0$.

###  II.2 Modified immittance function

Recall that the goal is to derive an expression for the immittance function of
the electrochemical system which necessitates an expression for the Laplace
transform of the current density in response to a voltage excitation across
the electrodes. In practice, an electrical stimulus (voltage or current) is
applied and the response (current or voltage) is measured for a discrete set
of frequencies, from which several information about the system can be
determined Barsoukov and Macdonald (2005). The system is characterized by its
transfer function $H(s)={\mathcal{L}{(v(t))}}/{\mathcal{L}(i(t))}$ or its
reciprocal. When a steady-state sinusoidal excitation is applied, and if the
system under test is also linear and the conditions of causality and time-
invariance are obeyed (Kramers–Kronig relationships), the transfer function
may be identified as the impedance $Z(j\omega)$ of the system, and its
reciprocal is the admittance function $Y(j\omega)$.

For analytical convenience, we consider here the simple case of ideally
polarized electrode submerged in dilute symmetric electrolyte with the ionic
charge of positive and negative ions being $|z_{{}_{+}}|=|z_{{}_{-}}|=z$ and
their diffusion coefficients being
$D_{{+_{\alpha}}}=D_{-_{\alpha}}=D_{\alpha}$. We also restrict our analysis to
one-dimensional case. In this case, the density of charge for any time $t$ and
location ${x}$ is $\rho({x},t)=zF(c_{{}_{+}}({x},t)-c_{{}_{-}}({x},t))$. Also,
in view of small perturbations of potential and charge, the sum of
compositions $c_{{}_{+}}({x},t)+c_{{}_{-}}({x},t)$ can be assumed to be
constant and always equal to $2c_{0}$. Thus, the current density at any cross
section of the system and at any time is obtained from Eq. 15 (or equivalently
Eq. 11 for the classical case) as Buck (1969); Bazant, Thornton, and Ajdari
(2004):

$i_{\alpha}(x,t)=-D_{\alpha}\frac{\partial\rho(x,t)}{\partial
x}-D_{\alpha}\kappa^{2}\epsilon\frac{\partial\phi(x,t)}{\partial x}$ (21)

where $\kappa$ is the reciprocal Debye thickness given by:

$\kappa=\sqrt{2c_{0}z^{2}F^{2}/(\epsilon RT)}$ (22)

Eq. 21 is the sum of two currents: a conduction current caused by
concentration gradients and a displacement current caused by the electric
field. Note that although we allow the potential and currents to be time-
varying, magnetic effects can be considered negligible in our situation
Bonnefont, Argoul, and Bazant (2001).

The domain’s boundary conditions we consider are such that the fluxes
${J}_{+_{\alpha}}={J}_{-_{\alpha}}=0$ at both electrodes for any time $t>0$
(i.e. no Faradaic current, perfect blocking electrodes), and
$\rho({x},0)=zF\left(c_{{}_{+}}({x},0)-c_{{}_{-}}({x},0)\right)=\rho_{0}$ and
$\rho({x},\infty)=\rho_{\infty}$ Buck (1969).

To solve this problem we proceed as follows. We assume for simplicity the same
values of delays $\tau_{1}$ and $\tau_{2}$, and the same order of fractional
derivative (i.e. $\alpha$) in Eq. 19 for describing the dynamics of both
positively and negatively-charged species. We write the fractional time
derivative of $\rho({x},t)=zF(c_{{}_{+}}({x},t)-c_{{}_{-}}({x},t))$ with the
use of Eqs. 19 and 15 (or equivalently the first-order time derivatives with
Eqs. 12 and 11 for the classical case) for both ionic species, which leads to
the following equation for the density of charge:

${}_{0}D_{t}^{\alpha}\rho(x,t+\tau_{1})=D_{\alpha}\frac{\partial^{2}\rho(x,t+\tau_{2})}{\partial
x^{2}}-D_{\alpha}\kappa^{2}\rho(x,t+\tau_{2})$ (23)

In detail this is done by first taking the divergence of the flux for both
species (Eq. 16 with separate delays), i.e.:

$\frac{\partial J_{+_{\alpha}}\left({x},t+\tau_{1}\right)}{\partial
x}=-D_{\alpha}\frac{\partial}{\partial
x}\left\\{\left[\frac{\partial}{\partial
x}+\frac{zF}{RT}\frac{\partial}{\partial
x}\phi\left({x},t+\tau_{2}\right)\right]c_{+}\left({x},t+\tau_{2}\right)\right\\}$
(24) $\frac{\partial{J}_{-_{\alpha}}\left({x},t+\tau_{1}\right)}{\partial
x}=-D_{\alpha}\frac{\partial}{\partial
x}\left\\{\left[\frac{\partial}{\partial
x}-\frac{zF}{RT}\frac{\partial}{\partial
x}\phi\left({x},t+\tau_{2}\right)\right]c_{-}\left({x},t+\tau_{2}\right)\right\\}$
(25)

Subtracting Eq. 25 from 24, multiplying by $zF$, and with the use of the
fractional continuity equation and
$2c_{0}=c_{+}\left({x},t\right)+c_{-}\left({x},t\right)$ we get:

${}_{0}D_{t}^{\alpha}\rho\left({x},t+\tau_{1}\right)=D_{\alpha}\left[\frac{\partial^{2}\rho\left({x},t+\tau_{2}\right)}{\partial
x^{2}}+\frac{z^{2}F^{2}}{RT}\frac{\partial}{\partial
x}\left[2c_{0}\frac{\partial\phi\left({x},t+\tau_{2}\right)}{\partial
x}\right]\right]$ (26)

With the Poisson equation, we finally obtain the result given in Eq. 23. Eq.
23 is similar to the classical PNP system Buck (1969):

$\frac{\partial\rho(x,t)}{\partial t}=D\frac{\partial^{2}\rho(x,t)}{\partial
x^{2}}-D\kappa^{2}\rho(x,t)$ (27)

when $\tau_{1}=\tau_{2}=0$. It can also be rewritten in dimensionless form as:

${}_{0}D_{\theta}^{\alpha}\rho(z,\theta+\theta_{1})=\frac{\partial^{2}\rho(z,\theta+\theta_{2})}{\partial
z^{2}}-\rho(z,\theta+\theta_{2})$ (28)

where $\theta=D_{\alpha}\kappa^{2}t^{1/\alpha}$,
$\theta_{j}=D_{\alpha}\kappa^{2}\tau_{j}^{1/\alpha}$, and $z=\kappa x$.

Now we apply the Laplace transform (defined as
$\tilde{f}(x,s)=\int_{0}^{\infty}f(x,t)e^{-st}\mathrm{d}t$) to both sides of
Eq. 23, considering zero initial conditions, to obtain:

$e^{-s\Delta\tau}s^{\alpha}\tilde{\rho}(x,s)=D_{\alpha}\frac{\mathrm{d}^{2}\tilde{\rho}(x,s)}{\mathrm{d}x^{2}}-D_{\alpha}\kappa^{2}\tilde{\rho}(x,s)$
(29)

Rearranging gives:

$\frac{\mathrm{d}^{2}\tilde{\rho}(x,s)}{\mathrm{d}x^{2}}-\tilde{\eta}_{\alpha}^{2}\tilde{\rho}(x,s)=0$
(30)

where

$\tilde{\eta}_{\alpha}^{2}=\frac{e^{-s\Delta\tau}s^{\alpha}}{D_{\alpha}}+\kappa^{2}$
(31)

It is straightforward to realize that when $\alpha=1$ and
$\Delta\tau=\tau_{2}-\tau_{1}=0$, Eq. 30 reduces to the equivalent classical
expression
${\mathrm{d}^{2}\tilde{\rho}}/{\mathrm{d}x^{2}}-\tilde{\eta}^{2}\tilde{\rho}=0$
where $\tilde{\eta}^{2}={s}/{D}+\kappa^{2}$ Gomez-Zamudio  _et al._ (2016).
The general solution of Eq. 30 is of the form:
$\tilde{\rho}(x,s)=\tilde{C_{1}}e^{-\tilde{\eta}_{\alpha}x}+\tilde{C_{2}}e^{\tilde{\eta}_{\alpha}x}$,
which with the conditions $\rho({x},0)=0$, $\rho(\infty,t)=0$ and
$\rho(0,t)=-\kappa^{2}\epsilon\phi_{0}(t)$ as used in Gomez-Zamudio et al.
Gomez-Zamudio  _et al._ (2016), we find:

$\tilde{\rho}(x,s)=-\kappa^{2}\epsilon\tilde{\phi}_{0}e^{-\tilde{\eta}_{\alpha}x}$
(32)

Now the Laplace transform applied to the Poisson equation (Eq. 20) with the
result of Eq 32 gives the following equation:

$\frac{\mathrm{d}^{2}\tilde{\phi}(x,s)}{\mathrm{d}x^{2}}=\kappa^{2}\tilde{\phi}_{0}e^{-\tilde{\eta}_{\alpha}x}$
(33)

which, after integrating once and using the condition that $\tilde{i}(0,s)=0$,
gives:

${\frac{\mathrm{d}\tilde{\phi}(x,s)}{\mathrm{d}x}=-\frac{\kappa^{2}\tilde{\phi}_{0}e^{-\tilde{\eta}_{\alpha}x}}{\tilde{\eta}_{\alpha}}+\frac{\tilde{\phi}_{0}(\kappa^{2}-\tilde{\eta}_{\alpha}^{2})}{\tilde{\eta}_{\alpha}}}$
(34)

Using Eq. 34 and the derivative of $\tilde{\rho}(x,s)$ (Eq. 32) w.r.t. $x$ we
obtain the following expression for the Laplace transform of the current:

$\tilde{i}_{\alpha}(x,s)=\frac{D_{\alpha}\kappa^{2}\epsilon\tilde{\phi}_{0}}{\tilde{\eta}_{\alpha}}\left(\kappa^{2}-\tilde{\eta}_{\alpha}^{2}\right)\left(e^{-\tilde{\eta}_{\alpha}x}-1\right)$
(35)

and thus a modified admittance function as:

$\displaystyle Y_{\alpha}(x,s)$
$\displaystyle=Z_{\alpha}^{-1}(x,s)=\frac{\tilde{i}_{\alpha}(x,s)}{\tilde{\phi}_{0}}$
$\displaystyle=\frac{D_{\alpha}\kappa^{2}\epsilon}{\tilde{\eta}_{\alpha}}\left(\kappa^{2}-\tilde{\eta}_{\alpha}^{2}\right)\left(e^{-\tilde{\eta}_{\alpha}x}-1\right)$
(36)

The expression of $Y_{\alpha}(x,s)$ incorporates both the fractional order
behavior of CPEs and a delay $\Delta\tau$ that limits the system bandwidth at
which energy storage takes place, as shall be shown in the next section. In
principle, the measurement of admittance from voltage excitation and current
response takes place at the boundaries of the EDLC system, and thus
$Y_{\alpha}(x,s)=Y_{\alpha}(x=x_{h},s)$ where $x_{h}$ denotes for instance the
thickness of the Helmholtz layer. This makes the admittance to be solely
dependent on the frequency in this case.

To find the real and imaginary part of $Y_{\alpha}(x,s)$, we write
$\tilde{\eta}_{\alpha}^{2}={e^{-s\Delta\tau}s^{\alpha}}/{D_{\alpha}}+\kappa^{2}=a+jb={(p+jq)^{2}}$
where:

$\displaystyle a$
$\displaystyle=\frac{\omega^{\alpha}\cos(\alpha\pi/2-\omega\Delta\tau)}{D_{\alpha}}+\kappa^{2}$
(37) $\displaystyle b$
$\displaystyle=\frac{\omega^{\alpha}\sin(\alpha\pi/2-\omega\Delta\tau)}{D_{\alpha}}$
(38)

which gives:

$\displaystyle p$
$\displaystyle=\frac{1}{\sqrt{2}}\sqrt{\sqrt{a^{2}+b^{2}}+a}$ (39)
$\displaystyle q$
$\displaystyle=\frac{\mathrm{sgn}\,{b}}{\sqrt{2}}\sqrt{\sqrt{a^{2}+b^{2}}-a}$
(40)

where $\mathrm{sgn}\,{b}=b/|b|$. Thus:

$\displaystyle\frac{\mathrm{Re}(Y_{\alpha})}{D\kappa^{2}\epsilon/(p^{2}+q^{2})}$
$\displaystyle=p(\kappa^{2}-p^{2}-q^{2})(e^{-px}\cos(qx)-1)$
$\displaystyle-q(\kappa^{2}+p^{2}+q^{2})e^{-px}\sin(qx)$ (41)
$\displaystyle\frac{\mathrm{Im}(Y_{\alpha})}{D\kappa^{2}\epsilon/(p^{2}+q^{2})}$
$\displaystyle=-q(\kappa^{2}+p^{2}+q^{2})(e^{-px}\cos(qx)-1)$
$\displaystyle-p(\kappa^{2}-p^{2}-q^{2})e^{-px}\sin(qx)$ (42)

from which the admittance phase angle is:

$\phi(Y_{\alpha})=\tan^{-1}[\mathrm{Im}(Y_{\alpha})/\mathrm{Re}(Y_{\alpha})]$
(43)

We can also express the results in terms of impedance such that:

$\displaystyle{\mathrm{Re}(Z_{\alpha})}={{\mathrm{Re}(Y_{\alpha})}}/{|Y_{\alpha}|^{2}}$
(44)
$\displaystyle{\mathrm{Im}(Z_{\alpha})}=-{{\mathrm{Im}(Y_{\alpha})}}/{|Y_{\alpha}|^{2}}$
(45) $\displaystyle\phi(Z_{\alpha})=-\phi(Y_{\alpha})$ (46)

which will be the focus of the simulation results below.

Figure 2: Plots of Eq. 36 in Nyquist impedance representation and phase angle
vs. log of frequency for (a)-(b) different value of $\Delta\tau$ and constant
values of $\alpha=0.5$, and (c)-(d) different values of $\alpha$ and constant
value of $\Delta\tau=0.01$ s. The parameters are $x=\text{1\,nm}$;
$T=\text{300\,K}$; $D=2\times 10^{-9}\,\text{m}^{2}\,\text{s}^{-1}$;
$z=\text{1}$; $c_{0}=\text{1\,mmol\,L}^{-1}$; $\epsilon=6.95\times 10^{-10}$ F
m-1, and the frequency points ($f=\omega/(2\pi)$) are logarithmically spaced
from $10^{-2}$ to $10^{2}$ Hz with 20 pts per decade

##  III Results and Discussion

Plots of parametric imaginary vs. real parts of impedance and impedance phase
angle vs. frequency using Eq. 36 are given in Fig. 2. The plots show the
effects of $\Delta\tau$ varied between 0 and 0.1 s (Figs. 2(a) and 2(b)) and
the effect of $\alpha$ varied between 1.0 and 0.3, i.e. subdiffusive transport
(Figs. 2(c) and 2(d)). We used the same parameters values as in Gomez-Zamudio
et al. Gomez-Zamudio  _et al._ (2016), i.e. $T=\text{300\,K}$;
$D_{\alpha}=2\times 10^{-9}\,\text{m}^{2}\,\text{s}^{-\alpha}$; $z=\text{1}$;
$c_{0}=\text{1\,mmol\,L}^{-1}$; $\epsilon=6.95\times 10^{-10}$ F m-1, and
frequency points logarithmically spaced from $10^{-2}$ to $10^{2}$ Hz (20 pts
per decade). As for $x$ we set it equal to 1 nm which is roughly the thickness
of the Helmholtz layer, and at which the current and thus the impedance are
computed.

As expected for the case of an ideal capacitor (i.e. $\Delta\tau=0$ and
$\alpha=1$), the Nyquist plot is vertical with the high-frequency intercept on
the real axis at practically zero Ohm. The impedance phase angle is constant
right at -90 deg. for all frequencies. Such results are shown in all
subfigures of Fig. 2 as a comparative reference.

Now as we vary the value of the time lag $\Delta\tau$ for a constant value of
$\alpha=0.5$ (Warburg regime), it is clear from Fig.2(b) that there is a
critical frequency above which the impedance phase angle increases
exponentially towards larger values than -45 deg., observed at low-
frequencies. This is controlled by the term $e^{-s\Delta\tau}$ in
$\tilde{\eta}_{\alpha}$ that appears in the expression of the modified
impedance function, which, in the time domain differential equations is
associated with $\Delta\tau=\tau_{2}-\tau_{1}$. It is important to mention
that the values of $\Delta\tau$ are positive in our case. This is different
from standard Catteneo’s formalism (purely due to diffusive transport, i.e.
$\tau_{2}=0$) because of the fact that the current in the PNP system is the
sum of both diffusion currents due to concentration gradients and drift
currents due to potential gradient. As a consequence, the net value of
$\Delta\tau$ can be thought of as the resulting time shift due to both
diffusion and drift. For the diffusion transport, it is the variation of
concentration that leads to the diffusion current after a certain time delay,
whereas the drift currents, from both positive and negative charges because of
the electric field, lead to changes in the concentration of the different
species after another time delay. It is clear that with the positive values of
$\Delta\tau$, the results are very comparable to real resistive-capacitive
devices as shown in Fig. 1 in which (i) both real and imaginary parts of the
impedance are frequency-dependent, but (ii) also is the impedance phase angle
especially at intermediate frequencies away from dc. This is not case for a
single CPE that is by definition known to exhibit a constant impedance phase
angle at all frequencies.

As for the effect of the fractional coefficient $\alpha$ for a fixed value of
$\Delta\tau$ (taken as 0.01 s) on the spectral impedance, Figs.2(c) and 2(d)
show clear angular deviation of the Nyquist plot and phase angle plot from
that of ideal capacitors as the values of $\alpha$ are decreased further away
from one. At the low frequencies limits, $\phi(Z)\to-\alpha\pi/2$ as expected
for a CPE. With $\alpha=1$, one recovers the normal electrodiffusion model and
thus the pure capacitive behavior. Again the presence of non-zero value for
$\Delta\tau$ in the expression of the modified impedance makes the system veer
towards a resistor as the frequency is increased.

In summary, the anomalous electrodiffusion problem described with fractional-
order calculus is known to give rise to the CPE behavior, which appears as a
straight line in the Nyquist impedance representation, or a constant impedance
phase for all frequencies. However, this is not the case for experimental
results of EDLCs for instance, which show a CPE behavior only for a limited
bandwidth at close-to-dc frequencies and a tendency towards resistive behavior
as the frequency is increased (Fig. 1). Lenzi et al. have just recently
demonstrated that with the use of fractional differential operators with non-
singular kernels (i.e., Atangana-Baleanu Atangana and Baleanu (2016)) it is
possible the retrieve the impedance response of a CPE in association with a
resistor, which is in very good agreement with the experiment. Nevertheless,
this approach relies on properly choosing a mathematical expression for the
kernel used in the fractional integro-differential equation, such as the
generalized two-parameter Mittag-Leffler function for instance, with which the
physical traceability may be blurred. Our modified PNP scheme with a time
shift between cause and effect, and with mono-order fractional time
derivatives also successfully captured such a resistive-capacitive behavior in
a compact and concise way. While the physical meaning of the fractional
exponent is still under debate, the use of delayed equations is believed to be
more advantageous to explain the physical mechanism in fractional EDLCs
Allagui  _et al._ (2020a).

##  IV Conclusion

In the standard electrokinetic PNP model, the Poisson equation relates the
free charge density to the Laplacian of the electric potential, and the
transport of dilute ionic species is governed by the fluxes of
electromigration and diffusion. The continuity equation connects the first-
order time derivative of ion concentrations to the divergence of the the total
flux. Here we proposed an extended fractional time derivative of time-shifted
concentrations of electrodiffusing species as a modified model for the PNP
system. The generalization of the integer-order time derivative to fractional-
order time derivative is commonly employed for describing the behavior of
complex system when integer-order models fail. Considering blocking electrodes
conditions, the model provided close-to-realistic description of typical EDLC
impedance profile at low-to-intermediate frequencies. With fractional
exponents different from one, the impedance phase angle at low frequencies is
shown to deviate from the response of ideal capacitors (i.e. CPE impedance),
whereas non-zero values of the time shift $\Delta\tau$ allowed the modeling of
the increased resistive behavior as the frequency is increased. The time shift
is due to the resulting net effect of both diffusion and drift time shifts
with respect to changes in concentration of diffusing species. In principle,
other boundary conditions such as Faradaic currents at the electrodes can be
investigated for simulating the impedance profiles of electrochemical sensors
or batteries.

## Acknowledgement

This work was supported by the NSF project #2126190 (C.W & A.A.)

## Data Availability

The data that supports the findings of this study are available within the
article.

## References

  * Allagui  _et al._ (2020a) A. Allagui, H. Alnaqbi, A. S. Elwakil, Z. Said, A. Hachicha, C. Wang, and M. A. Abdelkareem, “Fractional-order electric double-layer capacitors with tunable low-frequency impedance phase angle and energy storage capabilities,” Appl. Phys. Lett. 116, 013902 (2020a).
  * Allagui  _et al._ (2016a) A. Allagui, T. J. Freeborn, A. S. Elwakil, and B. J. Maundy, “Reevaluation of performance of electric double-layer capacitors from constant-current charge/discharge and cyclic voltammetry,” Sci. Rep. 6 (2016a).
  * Allagui  _et al._ (2018a) A. Allagui, T. J. Freeborn, A. S. Elwakil, M. E. Fouda, B. J. Maundy, A. G. Radwanh, Z. Said, and M. A. Abdelkareem, “Review of fractional-order electrical characterization of supercapacitors,” J. Power Sources 400 (2018a).
  * Allagui  _et al._ (2016b) A. Allagui, A. S. Elwakil, B. J. Maundy, and T. J. Freeborn, “Spectral capacitance of series and parallel combinations of supercapacitors,” ChemElectroChem 3, 1429–1436 (2016b).
  * Allagui, Elwakil, and Psychalinos (2021) A. Allagui, A. S. Elwakil, and C. Psychalinos, “Decoupling the magnitude and phase in a constant phase element,” J. Electroanal. Chem. 888, 115153 (2021).
  * Bisquert  _et al._ (1998) J. Bisquert, G. Garcia-Belmonte, P. Bueno, E. Longo, and L. Bulhoes, “Impedance of constant phase element (cpe)-blocked diffusion in film electrodes,” J. Electroanal. Chem. 452, 229–234 (1998).
  * Allagui, Elwakil, and Fouda (2021) A. Allagui, A. S. Elwakil, and M. E. Fouda, “Revisiting the time-domain and frequency-domain definitions of capacitance,” IEEE Trans. Electron Devices 68 (2021).
  * Allagui and Fouda (2021) A. Allagui and M. E. Fouda, “Inverse problem of reconstructing the capacitance of electric double-layer capacitors,” Electrochim. Acta , 138848 (2021).
  * Allagui  _et al._ (2020b) A. Allagui, D. Zhang, I. Khakpour, A. S. Elwakil, and C. Wang, “Quantification of memory in fractional-order capacitors,” J. Phys. D 53 (2020b).
  * Allagui, Zhang, and Elwakil (2018) A. Allagui, D. Zhang, and A. S. Elwakil, “Short-term memory in electric double-layer capacitors,” Appl. Phys. Lett. 113, 253901–5 (2018).
  * Baboukani  _et al._ (2019) A. R. Baboukani, I. Khakpour, V. Drozd, A. Allagui, and C. Wang, “Single-step exfoliation of black phosphorus and deposition of phosphorene via bipolar electrochemistry for capacitive energy storage application,” J. Mater. Chem. A 7 (2019).
  * Khakpour  _et al._ (2019) I. Khakpour, A. R. Baboukani, A. Allagui, and C. Wang, “Bipolar exfoliation and in-situ deposition of high-quality graphene for supercapacitor application,” ACS Appl. Energy Mater. 2, 4813–4820 (2019).
  * Allagui  _et al._ (2021) A. Allagui, A. R. Baboukani, A. S. Elwakil, and C. Wang, “Electrochemical stability analysis of red phosphorus-based anode for lithium-ion batteries,” Electrochim. Acta , 139149 (2021).
  * Cheng  _et al._ (2018) C.-S. Cheng, H. S.-H. Chung, R. W.-H. Lau, and K. Y.-W. Hong, “Time-domain modeling of constant phase elements for simulation of lithium battery behavior,” IEEE Trans. Power Electron. 34, 7573–7587 (2018).
  * Xiong  _et al._ (2018) R. Xiong, J. Tian, W. Shen, and F. Sun, “A novel fractional order model for state of charge estimation in lithium ion batteries,” IEEE Trans. Veh. Technol. 68, 4130–4139 (2018).
  * Cuervo-Reyes  _et al._ (2015) E. Cuervo-Reyes, C. P. Scheller, M. Held, and U. Sennhauser, “A unifying view of the constant-phase-element and its role as an aging indicator for li-ion batteries,” J. Electrochem. Soc. 162, A1585 (2015).
  * Zhang  _et al._ (2019) D. Zhang, A. Allagui, A. S. Elwakil, A. M. Nassef, H. Rezk, J. Chengi, and W. C. Choy, “On the modeling of dispersive transient photocurrent response of organic solar cells,” Org. Electron. 70, 42–47 (2019).
  * Zhang  _et al._ (2020) D. Zhang, A. Allagui, A. S. Elwakil, Z. Yan, and H. Lu, “Active circuit model of low-frequency behavior in perovskite solar cells,” Org. Electron. 85, 105804 (2020).
  * Borini  _et al._ (2013) S. Borini, R. White, D. Wei, M. Astley, S. Haque, E. Spigone, N. Harris, J. Kivioja, and T. Ryhanen, “Ultrafast graphene oxide humidity sensors,” ACS nano 7, 11166–11173 (2013).
  * Song  _et al._ (2020) J. Song, Y. Li, F. Yin, Z. Zhang, D. Ke, D. Wang, Q. Yuan, and X.-E. Zhang, “Enhanced electrochemical impedance spectroscopy analysis of microbial biofilms on an electrochemically in situ generated graphene interface,” ACS sensors 5, 1795–1803 (2020).
  * Meunier  _et al._ (2020) C. J. Meunier, J. D. Denison, G. S. McCarty, and L. A. Sombers, “Interpreting dynamic interfacial changes at carbon fiber microelectrodes using electrochemical impedance spectroscopy,” Langmuir 36, 4214–4223 (2020).
  * Allagui and Elwakil (2021) A. Allagui and A. S. Elwakil, “Possibility of information encoding/decoding using the memory effect in fractional-order capacitive devices,” Sci. Rep. 11, 1–7 (2021).
  * Podlubny (1998) I. Podlubny,  _Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications_ (Elsevier, 1998).
  * Fouda  _et al._ (2016) M. E. Fouda, A. S. Elwakil, A. G. Radwan, and A. Allagui, “Power and energy analysis of fractional-order electrical energy storage devices,” Energy 111, 785–792 (2016).
  * Allagui  _et al._ (2018b) A. Allagui, A. S. Elwakil, Z. Said, M. A. Abdelkareem, and D. Zhang, “Band-pass filter and relaxation oscillator using electric double-layer capacitor,” ChemElectroChem 5 (2018b).
  * Allagui  _et al._ (2017) A. Allagui, Z. Said, M. A. Abdelkareem, A. S. Elwakil, M. Yang, and H. Alawadhi, “DC and AC performance of graphite films supercapacitors prepared by contact glow discharge electrolysis,” J. Electrochem. Soc. 164, A2539–A2546 (2017).
  * Kosztołowicz and Lewandowska (2009) T. Kosztołowicz and K. D. Lewandowska, “Hyperbolic subdiffusive impedance,” J. Phys. A: Math. Theor. 42, 055004 (2009).
  * Kant, Kumar, and Yadav (2008) R. Kant, R. Kumar, and V. K. Yadav, “Theory of anomalous diffusion impedance of realistic fractal electrode,” J. Phys. Chem. C 112, 4019–4023 (2008).
  * Maritan, Stella, and Toigo (1989) A. Maritan, A. Stella, and F. Toigo, “Anomalous warburg impedance and universal surface magnetic exponent for gaussian models in the presence of fractal boundaries,” Phys. Rev. B 40, 9269 (1989).
  * Zhang  _et al._ (2015) Y. Zhang, Y. Chen, M. Yan, and F. Chen, “Reconstruction of relaxation time distribution from linear electrochemical impedance spectroscopy,” J. Power Sources 283, 464–477 (2015).
  * Barbero (2017) G. Barbero, “Theoretical interpretation of warburg’s impedance in unsupported electrolytic cells,” Physical Chemistry Chemical Physics 19, 32575–32579 (2017).
  * Kilic, Bazant, and Ajdari (2007) M. S. Kilic, M. Z. Bazant, and A. Ajdari, “Steric effects in the dynamics of electrolytes at large applied voltages. ii. modified poisson-nernst-planck equations,” Phys. Rev. E 75, 021503 (2007).
  * Schmuck and Bazant (2015) M. Schmuck and M. Z. Bazant, “Homogenization of the poisson–nernst–planck equations for ion transport in charged porous media,” SIAM J. Appl. Math. 75, 1369–1401 (2015).
  * Alexe-Ionescu  _et al._ (2020) A. Alexe-Ionescu, G. Barbero, L. Evangelista, and E. Lenzi, “Current–voltage characteristics and impedance spectroscopy: Surface conduction and adsorption–desorption effects in electrolytic cells,” J. Phys. Chem. C 124, 3150–3158 (2020).
  * Koklu, Ajaev, and Beskok (2019) A. Koklu, V. Ajaev, and A. Beskok, “Self-similar response of electrode polarization for binary electrolytes in parallel plate capacitor systems,” Analytical chemistry 91, 11231–11239 (2019).
  * Compte and Metzler (1997) A. Compte and R. Metzler, “The generalized cattaneo equation for the description of anomalous transport processes,” J. Phys. A: Math. Gen. 30, 7277 (1997).
  * Lewandowska and Kosztołowicz (2008) K. D. Lewandowska and T. Kosztołowicz, “Application of generalized cattaneo equation to model subdiffusion impedance,” Acta Phys. Pol. B 39 (2008).
  * Ruggeri, Muracchini, and Seccia (1990) T. Ruggeri, A. Muracchini, and L. Seccia, “Shock waves and second sound in a rigid heat conductor: A critical temperature for naf and bi,” Physical review letters 64, 2640 (1990).
  * Joseph and Preziosi (1989) D. D. Joseph and L. Preziosi, “Heat waves,” Rev. Mod. Phys. 61, 41–73 (1989).
  * Kosztołowicz (2014) T. Kosztołowicz, “Cattaneo-type subdiffusion-reaction equation,” Physical Review E 90, 042151 (2014).
  * Garcia-Colin and Olivares-Robles (1995) L. Garcia-Colin and M. A. Olivares-Robles, “Hyperbolic type transport equations,” Physica A: Statistical Mechanics and its Applications 220, 165–172 (1995).
  * Qi and Jiang (2011) H. Qi and X. Jiang, “Solutions of the space-time fractional cattaneo diffusion equation,” Physica A: Statistical Mechanics and its Applications 390, 1876–1883 (2011).
  * Gomez-Zamudio  _et al._ (2016) F. Gomez-Zamudio, R. Antano-Lopez, A. Rodriguez-Lopez, T. Perez, and E. R. Larios-Duran, “Development of the electrochemical impedance response of ideally polarized interfaces based on transport phenomena laws through the nernst-planck-poisson equation linearized by the debye-falkenhagen approximation,” Int. J. Electrochem. Sci. 11, 475 – 485 (2016).
  * Lenzi  _et al._ (2021) E. Lenzi, L. Guilherme, B. da Silva, A. Koltun, L. Evangelista, and R. Zola, “Anomalous diffusion and electrical impedance response: Fractional operators with singular and non-singular kernels,” Commun. Nonlinear Sci. Numer. Simul. 102, 105907 (2021).
  * Hilfer (1995) R. Hilfer, “Foundations of fractional dynamics,” Fractals 3, 549–556 (1995).
  * Manzanares and Kontturi (2003) J. A. Manzanares and K. Kontturi, “Diffusion and migration,” Encyclopedia of electrochemistry 2, 81–121 (2003).
  * Hilfer (2000) R. Hilfer, “Fractional diffusion based on riemann-liouville fractional derivatives,” The Journal of Physical Chemistry B 104, 3914–3917 (2000).
  * Tzou (1995) D. Y. Tzou, “A unified field approach for heat conduction from macro-to micro-scales,” J. Heat Transfer 117, 8–16 (1995).
  * Górska (2021) K. Górska, “Integral decomposition for the solutions of the generalized cattaneo equation,” Physical Review E 104, 024113 (2021).
  * Lenzi, Evangelista, and Barbero (2009) E. Lenzi, L. Evangelista, and G. Barbero, “Fractional diffusion equation and impedance spectroscopy of electrolytic cells,” The Journal of Physical Chemistry B 113, 11371–11374 (2009).
  * Kosztołowicz (2019) T. Kosztołowicz, “Model of anomalous diffusion-absorption process in a system consisting of two different media separated by a thin membrane,” Physical Review E 99, 022127 (2019).
  * Barsoukov and Macdonald (2005) E. Barsoukov and J. R. Macdonald,  _Impedance Spectroscopy Theory, Experiment, and Applications_ (Wiley Online Library, 2005).
  * Buck (1969) R. Buck, “Diffuse layer charge relaxation at the ideally polarized electrode,” J Electroanal Chem Interfacial Electrochem. 23, 219–240 (1969).
  * Bazant, Thornton, and Ajdari (2004) M. Z. Bazant, K. Thornton, and A. Ajdari, “Diffuse-charge dynamics in electrochemical systems,” Physical review E 70, 021506 (2004).
  * Bonnefont, Argoul, and Bazant (2001) A. Bonnefont, F. Argoul, and M. Z. Bazant, “Analysis of diffuse-layer effects on time-dependent interfacial kinetics,” J. Electroanal. Chem. 500, 52–61 (2001).
  * Atangana and Baleanu (2016) A. Atangana and D. Baleanu, “New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model,” arXiv preprint arXiv:1602.03408 (2016).