# Electrochemical impedance spectroscopy beyond linearity and stationarity —
a critical review

Noël Hallemans<EMAIL_ADDRESS>David Howey  Alberto Battistel  Nessa
Fereshteh Saniee  Federico Scarpioni  Benny Wouters  Fabio La Mantia  Annick
Hubin  Widanalage Dhammika Widanage  John Lataire  Research Group Fundamental
Electricity and Instrumentation, Vrije Universiteit Brussel, Pleinlaan 2, 1050
Brussels, Belgium  WMG, University of Warwick, Coventry, 7AL CV4, UK  Battery
Intelligence Lab, Department of Engineering, University of Oxford, OX1 3PJ, UK
The Faraday Institution, Quad One, Harwell Science and Innovation Campus,
Didcot, UK  Institute of Technical Medicine, Furtwangen University, Jakob-
Kienzle-Strasse 17, Villingen-Schwenningen 78054, Germany  Fraunhofer
Institute for Manufacturing Technology and Advanced Materials IFAM, Wiener
Strasse 12, Bremen 28359, Germany  Research Group Electrochemical and Surface
Engineering, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium
Energy Storage and Energy Conversion Systems, Bremen University, Wiener
Strasse 12, Bremen 28359, Germany

###### Abstract

Electrochemical impedance spectroscopy (EIS) is a widely used experimental
technique for characterising materials and electrode reactions by observing
their frequency-dependent impedance. Classical EIS measurements require the
electrochemical process to behave as a linear time-invariant system. However,
electrochemical processes do not naturally satisfy this assumption: the
relation between voltage and current is inherently nonlinear and evolves over
time. Examples include the corrosion of metal substrates and the cycling of
Li-ion batteries. As such, classical EIS only offers models linearised at
specific operating points. During the last decade, solutions were developed
for estimating nonlinear and time-varying impedances, contributing to more
general models. In this paper, we review the concept of impedance beyond
linearity and stationarity, and detail different methods to estimate this from
measured current and voltage data, with emphasis on frequency domain
approaches using multisine excitation. In addition to a mathematical
discussion, we measure and provide examples demonstrating impedance estimation
for a Li-ion battery, beyond linearity and stationarity, both while resting
and while charging.

###### keywords:

EIS , dynamic EIS , NLEIS , impedance , multisine , nonlinearity ,
nonstationarity , frequency domain , Li-ion , battery

††journal: Electrochimica Acta

Electrochemistry studies processes at electrode/electrolyte interfaces. These
processes involve the movement of charged species (ions or electrons),
generating a current flowing through a cell and a voltage drop over its
electrodes. Diverse noninvasive techniques relying on current and voltage
measurements have been developed for studying these processes. Typically, one
of the two quantities is kept constant, swept, or oscillated, while the other
quantity’s response is recorded. Some of the most widely used techniques are
linear sweep voltammetry, constant current chrono-potentiometry, constant-
potential chronoamperometry and electrochemical impedance spectroscopy (EIS).
In EIS [1, 2, 3, 4, 5], the dynamics of electrochemical processes are studied
by means of the impedance response to the applied current or voltage, and the
term ‘spectroscopy’ refers to the frequency dependency. It is worth
emphasising that EIS is a nonparametric data-driven technique, i.e. the
impedance is solely computed relying on current and voltage data, without
prior knowledge about the governing equations as in physics-based modelling
[6, 7, 8].

During the last two decades, EIS has gained wide popularity thanks to its
accessible implementation and broad applicability to, among others, corrosion
[9, 10, 11], batteries [12, 13, 14, 15, 16, 17, 18, 19, 20, 21], and fuel
cells [22, 23]. Nowadays, EIS is available in many commercial cyclers and
potentiostats, where a user decides on a set of frequencies and the device
measures the complex impedance values at these frequencies.

From a system theoretical perspective, the impedance is a special case of a
transfer function, which is a model for a linear time-invariant (LTI)
dynamical system [24]. In system theory, dynamical systems are systems with
memory, that is, systems defined through differential equations or,
equivalently, convolution operators. This is opposed to static systems, where
the output is simply a static function of the input. This interpretation of
the term dynamic should not be confused with its use in electrochemistry to
denote time-variation. In this article, we use the system theory convention.
Models relate the output of the system to its input, which, in the EIS case,
is the current through and the voltage over the electrodes. In _galvanostatic_
experiments the current is the input and voltage the output, while in
_potentiostatic_ experiments it is the other way around. In the remainder of
this text, EIS experiments satisfying the assumptions of LTI systems will be
referred to as ‘classical EIS’. Estimating transfer functions from input and
output data of LTI systems is a thoroughly studied problem in the field of
system identification [25, 26].

However, for electrochemical systems it is well known that (i) the relation
between current and voltage is generally nonlinear, e.g. expressed by
exponential functions in the Butler-Volmer equation, and (ii) the behaviour of
the process may evolve over time. Li-ion batteries, for example, provide such
time-variation on different time-scales. On a large time-scale, the impedance
of a fresh and an aged cell are different [12, 18]. On a shorter time-scale,
the impedance of a fully charged and a fully discharged cell are also
different [13]. Accordingly, with classical EIS, one represents
electrochemical processes which are by nature nonlinear and nonstationary with
a model for linear and stationary systems. While classical EIS may be a non-
ideal approximation for such cases, important information about the process is
nonetheless revealed by measuring within forced constraints of linearity and
stationarity at specific operating points. Linearity is achieved by applying a
small excitation amplitude such that the behaviour of the process is
linearised in a certain operating region. Stationarity is achieved by
measuring at a time, and over a timespan, when the process is and remains in
steady-state. However, these experimental conditions are very restrictive. How
can we obtain information about the nonlinear behaviour of a system when it is
only possible to perform experiments under linear constraints? How can we
study a battery _while_ it is charging or discharging? How can we study
anodising _while_ a protective layer is forming?

To relax these limitations, the concept of impedance needs extensions. The
required extensions to model systems beyond the linearity and stationarity
constraints have been studied in system theory. Nonlinear time-invariant
(NLTI) systems are commonly studied by _Volterra series_ [27, 28]. These are
convolution operators that are able to capture dynamical nonlinear behaviour,
where the transfer function is extended to so-called generalised transfer
functions. Similarly, nonlinear impedances are studied in nonlinear EIS
(NLEIS) [29, 30, 31, 32, 33, 34, 35]. Assuming that the nonlinear behaviour of
an electrochemical system can be captured by a Volterra series, the behaviour
of the process can be split into a linear part and purely nonlinear part [36].
The model for the linear part is denoted as the best linear approximation
(BLA) of the system. When the nonlinearities are small enough compared to the
linear behaviour, use of the BLA for describing the system is justified.

Nonstationary systems have been studied by extending the transfer function to
be a time-varying transfer function to describe a linear time-varying (LTV)
system [37]. The latter is a transfer function depending on both time and
frequency, expressing the evolution of the transfer function over time.
Similarly, the impedance can be extended to a _time-varying_ impedance. As an
analogy, the impedance can be considered to be like a photograph where as in
the first days of photography—think mid 19th century—the subject should remain
stationary during a certain exposure time for accumulating light on a sheet,
but the time-varying impedance can be seen as a movie of a subject during an
activity.

Over the years, techniques have been developed to _detect_ nonstationarities
and nonlinearities in measured data, where it was shown that frequency domain
identification techniques with multisine excitations are advantageous [38, 39,
40, 41, 42]. Recently, with increasing computation power, different techniques
have been developed and refined for unravelling time-varying impedance from
data [43, 44, 45, 46, 47, 48, 49]. Time-varying impedance data have already
successfully been obtained for a wide variety of electrochemical processes,
including organic coatings [50], electrorefining [51], hydrogen evolution
reactions [52], nickel hexacyanoferrate thin films [53], electrochemical
double layer capacitors [54], charging/discharging Li-ion batteries [55, 56,
57, 58, 59], and Li-plating in batteries [60, 61, 62].

Commonly, drift signals may appear (for instance the slow voltage increase
when charging a battery) during _in operando_ electrochemical measurements.
These drift signals prohibit measuring impedance data at low frequencies. A
method for removing drift signals has been developed [63], and successfully
applied to electrorefining [64], anodising [65], and Li-ion batteries [58,
59], such that the time variations of the low-frequency impedance can also be
estimated.

In this review article, we detail the required mathematical concepts
associated with impedance and their extensions beyond linearity and
stationarity. Key concepts are supported with illustrations obtained from
simulations and real-life measurements. Experiments are performed on a
pristine commercially available Samsung 48X Li-ion battery. This is a
$4.8\,\mathrm{Ah}$ $21$ $700$ cylindrical cell format with cathodes based on
lithiated metal oxide (Co, Ni, Al) and anodes based on graphite and blended
Si. The impedance is measured at different temperatures, while resting and
while charging. We opt for this case-study on Li-ion batteries since EIS is
becoming a popular tool for characterising batteries, diagnosing state-of-
health, and developing smart charging protocols. These compelling applications
are also discussed as a motivation to perform EIS beyond linearity and
stationarity. Moreover, the measurements are obtained using a commercial
potentiostat (Gamry Interface 5000E), showcasing the practical accessibility
of the discussed modelling techniques.

This article is structured as follows. First, we give a motivational example
on how impedance data beyond linearity and stationarity is promising for
battery aging diagnostics and smart charging protocols (Section 1). Next, we
define what we mean with a model for the electrochemical system (Section 2).
Then, we revisit classical EIS (Section 3), with emphasis on the limiting
constraints of linearity and stationarity. The choice between single-sine and
multisine excitations is discussed in depth. Then, we formally introduce
nonlinear and nonstationary models for electrochemical measurements through,
respectively, Volterra series and time-varying impedances (Section 4). The
Volterra series is linked to NLEIS and the BLA. Next, we detail the
experimental procedure in measuring current and voltage time-series for proper
impedance measurements (Section 5). The estimation of classical impedance data
in the frequency domain from periodic and random excitations is discussed in
Section 6. Then, we detail how the linearity and stationarity constraints can
be assessed, and nonlinearities and nonstationarities are detected by
observing the current and voltage spectra under odd random phase multisine
excitations (Section 7). Obtaining nonlinear impedance data and the BLA is
discussed in Section 8. The extraction of time-varying impedance data from the
collected current and voltage data is studied through different relevant
methods in Section 9. This has been studied in Szekeres et al. [66] also,
however, here a deeper mathematical foundation is given. In Section 10, the
performed illustrative experiments on Li-ion batteries are discussed. Finally,
conclusions are drawn and an outlook is given in Section 11.

##  1 Overview of applications of EIS for batteries

Before we get into the details of how impedance works, let us first motivate
the topic and reflect on why impedance is useful for solving some of
electrochemistry’s crucial research problems, and more importantly, why
measuring impedance beyond linearity and stationarity is promising. We do this
for the compelling case of Li-ion batteries. Here, some relevant research
problems include state-of-health (SOH) prognostics [67, 68, 69, 70, 71, 72,
73, 74, 75, 76] and smart charging [77, 78, 79, 80]. For SOH prognostics, EIS
is a powerful non-invasive tool [81]. It has been shown that classical
impedance data, mapped onto equivalent circuit model (ECM) parameters,
contains important information about degradation mechanisms [16, 73, 82].
Moreover, classical impedance data is an informative input to machine learning
algorithms to predict the remaining-useful-life (RUL) of batteries [19, 83].
Jones et al. [21, Table 1] demonstrate that the EIS data state representation
performs better than other state representations for SOH forecasting. Bizeray
et al. [84] also show that physical parameters can be identified from
classical impedance data, allowing simulation of the battery using physics-
based models such as the single-particle model (SPM). This is an non-invasive
way of parametrising cells, being more practical than tearing down cells [7].

Impedance data beyond linearity and stationarity has the potential to improve
battery SOH diagnostics since it contains additional information compared to
classical impedance data. Leveraging nonlinear and time-varying impedance data
has already been carried out for detecting Li-plating [85, 62] (an important
degradation mechanism in Li-ion batteries [70, 86]). Moreover, Kirk et al.
[34] demonstrate that nonlinear impedance data contributes to the
identifiability of physical SPM parameters.

Impedance data beyond linearity and stationarity also has the potential to
improve smart charging protocols. Katzer et al. [87] propose an adaptive fast
charging protocol relying on impedance-based detection of Li-plating [62]. In
Zhu et al. [58], we track the charge transfer resistance while charging, which
can be obtained from time-varying impedance data, and propose to adapt the
charge profile based on this time-varying charge transfer resistance.

A critical reader may argue that EIS experiments are expensive, and cannot
easily be implemented in battery management systems (BMS). However, different
solutions have been developed to implement low-cost measurement apparatus [88,
89, 90]. Here multisine excitations and frequency domain model estimation are
promising tools.

##  2 Modelling electrochemical systems

Electrochemical processes are often studied by modelling the relation between
the current $i(t)$ through, and the voltage $v(t)$ over the electrodes, where
$t$ denotes continuous time. However, this relation may also be dependent on
external parameters $p(t)$ such as the ambient temperature, the rotation rate
of the electrodes, the pressure, the concentration distribution on the
electrodes, etc.

Two types of experiments are common to measure the relation between current
and voltage. In _galvanostatic_ experiments, a current is applied and the
voltage is measured. In a general setting, the impedance is modelled as an
operator $\mathcal{G}\\{\cdot\\}$ acting on the current and external
parameters,

$\displaystyle v(t)=\mathcal{G}\\{i(t),p(t)\\}.$ (1)

In _potentiostatic_ experiments, the measurements are performed the other way
around, that is, through an operator $\mathcal{P}$,
$i(t)=\mathcal{P}\\{v(t),p(t)\\}$. Here, the notations for galvanostatic
experiments are chosen, that is, the current $i(t)$ is the excitation and the
voltage $v(t)$ is the response (appropriate for low impedance devices such as
batteries).

As mentioned in the introduction, the operators above are dynamical, that is,
operators with memory, or also called _convolution operators_ , not only
acting on the present time, but also making use of past information.

In what follows, it is detailed how the operator $\mathcal{G}$ can be modelled
by an impedance. We first consider the system under the constraints of
linearity and stationarity (Section 3), and proceed by introducing models
beyond these hard constraints (Section 4).

##  3 Classical EIS revisited

Figure 1: Equivalent electrical schematic of an electrochemical system under
LTI constraints.

###  3.1 The constraints of classical EIS

In classical EIS experiments, the external parameters $p(t)$ are assumed
constant during the experiment and the generic model (1) is simplified to

$\displaystyle v(t)=\mathrm{OCV}+\underbrace{Z\\{i(t)\\}}_{v_{Z}(t)},$ (2)

where OCV is the open circuit voltage, assumed constant, $Z$ is the classical
impedance operator and $v_{Z}(t)$ is the voltage over the impedance. An
equivalent electrical circuit representing this model is shown in Fig. 1. It
is assumed that operator $Z$ satisfies the constraints of LTI systems, that
is, linearity and stationarity, and also causality.

##### Linearity

The operator $Z$ is a linear operator, satisfying additivity and homogeneity,
respectively,

$\displaystyle Z\\{i_{1}(t)+i_{2}(t)\\}=Z\\{i_{1}(t)\\}+Z\\{i_{2}(t)\\}$ (3a)
$\displaystyle Z\\{\alpha i(t)\\}=\alpha
Z\\{i(t)\\}\quad\forall\alpha\in\mathbb{R}.$ (3b)

As such, when the current $i$ doubles, the voltage $v_{Z}$ over the impedance
also doubles.

##### Stationarity (time-invariance)

A stationary system is a system whose behaviour does not change when shifted
in time. Accordingly, the operator $Z$ is independent of the time at which the
excitation is applied:

$\displaystyle Z\\{i(t-\tau)\\}=v_{Z}(t-\tau)\quad\forall\tau\in\mathbb{R}.$
(4)

##### Causality

The response of the system is totally determined by the excitation. As a
consequence, the response to an excitation cannot precede the excitation.

For _potentiostatic_ experiments under LTI constraints, the excitation-
response relation yields,

$\displaystyle i(t)=Y\\{\underbrace{v(t)-\mathrm{OCV}}_{v_{Z}(t)}\\},$ (5)

where $Y$ is called the admittance operator, satisfying the same conditions as
the impedance operator $Z$. Notations in this paper can, hence, be converted
into potentiostatic experiments by swapping $i$ and $v_{Z}$, and replacing the
impedance $Z$ by admittance $Y$.

When charge transfer is the rate-determining step, the static relation between
current and voltage of electrochemical reactions is described by the Butler-
Volmer equation [91],

$\displaystyle
i=j_{0}S\left(\exp\left(\frac{\alpha_{a}nF}{R\mathrm{T}}v_{\mathrm{Z}}\right)-\exp\left(-\frac{\alpha_{c}nF}{R\mathrm{T}}v_{\mathrm{Z}}\right)\right),$
(6)

with $j_{0}$ the exchange current density, $S$ the surface area of the
electrode, $\mathrm{T}$ the absolute temperature in Kelvin, $n$ the number of
electrons, $F$ the Faraday constant, $R$ the universal gas constant,
$\alpha_{a}$ the anodic charge transfer coefficient, and $\alpha_{c}$ the
cathodic charge transfer coefficient. When assuming
$\alpha_{a}=\alpha_{c}=0.5$, this equation can be rewritten to obtain the
overpotential as a function of the current,

$\displaystyle
v_{\mathrm{Z}}=\frac{2R\mathrm{T}}{nF}\sinh^{-1}\left(\frac{i}{2j_{0}S}\right).$
(7)

The top left plot of Fig. 2 shows this Butler-Volmer relation between current
and voltage as a dotted line. This relation is obviously not linear. However,
the linearity constraint can approximately be satisfied by choosing the
magnitude of the excitation signal in a specific range. This can be seen by
expanding (7) as a Taylor series around $i=0$,

$\displaystyle v_{\mathrm{Z}}=\underbrace{\frac{\partial
v_{\mathrm{Z}}}{\partial i}\bigg{\rvert}_{i=0}i}_{\text{linear
term}}+\frac{\partial^{2}v_{\mathrm{Z}}}{\partial
i^{2}}\bigg{\rvert}_{i=0}i^{2}+\frac{\partial^{3}v_{\mathrm{Z}}}{\partial
i^{3}}\bigg{\rvert}_{i=0}i^{3}+\ldots$ (8)

When the current $i$ is small enough, the linear term will dominate the higher
order terms and linearity can be assumed. A rule of thumb for ensuring
linearity is that the voltage deviation should not be larger than
$15\text{\,}\mathrm{m}\mathrm{V}$ [92]. An illustration of the linearisation
of the Butler-Volmer equation is also shown in Fig. 2. A small amplitude
sinusoidal current excitation centered around zero (blue) is applied to the
Butler-Volmer equation (dotted black line), and the response is the voltage
centered around the value OCV (red). For an excitation with small amplitude,
the Butler-Volmer equation is quasi-linear in the excited range (black line).

The stationarity constraint, on the other hand, is satisfied by driving the
electrochemical system in steady-state, and applying a zero-mean current
excitation to remain in steady-state. For a battery, for instance, applying a
current excitation with a positive mean value would charge the battery and
cause nonstationary behaviour. It is also necessary that the external
parameters $p(t)$ remain constant during the experiment. A significant change
in ambient temperature, for instance, might jeopardise the stationarity
constraint.

Figure 2: Illustration of the linearisation of the Butler-Volmer equation (7)
for a small amplitude sinusoidal excitation.

The relation between the excitation and response of LTI systems is well
documented in system theory [24]. The response $v_{\mathrm{Z}}(t)$ of an LTI
system is commonly modelled by the _convolution_ of the impulse response
function $z(t)$ with the excitation $i(t)$,

$\displaystyle
v_{Z}(t)=\int_{-\infty}^{\infty}z(\tau)i(t-\tau)\mathrm{d}\tau.$ (9)

The impulse response function is the response to a Dirac pulse. Note that (9)
satisfies (3) and (4). Often, the fact that a convolution in the time domain
becomes a product in the frequency domain is exploited to rewrite (9) in a way
where the frequency dependent impedance $Z(\omega)$ appears,

$\displaystyle v_{Z}(t)=\mathcal{F}^{-1}\\{Z(\omega)I(\omega)\\}.$ (10)

Here, the impedance $Z(\omega)$ is defined as the Fourier transform of the
impulse response function $z(t)$, $\mathcal{F}^{-1}\\{\cdot\\}$ is the inverse
Fourier transform operator, and $I(\omega)$ the Fourier transform of the
current. The Fourier and inverse Fourier transforms are defined in A. Recall
that the angular frequency $\omega$ is related to the frequency $f$ as
$\omega=2\pi f$.

The voltage over the electrochemical system is then

$\displaystyle v(t)=\mathrm{OCV}+\mathcal{F}^{-1}\\{Z(\omega)I(\omega)\\}.$
(11)

Accordingly, the impedance is given by the _ratio of the Fourier transforms of
voltage and current_ ,

$\displaystyle Z(\omega)=\frac{V(\omega)}{I(\omega)}\qquad\omega\neq 0.$ (12)

Note that the impedance is not expressed at DC. The impedance at DC becomes
infinite in magnitude and with a purely imaginary phase because the linearised
OCV behaves like a capacitor.

As for any frequency response function, the impedance is a complex valued
function, often denoted as

$\displaystyle Z(\omega)=Z_{\mathrm{r}}(\omega)+jZ_{\mathrm{j}}(\omega),$ (13)

with $j$ the imaginary unit ($j^{2}=-1$), and $Z_{\mathrm{r}}(\omega)$ and
$Z_{\mathrm{j}}(\omega)$ the real and imaginary parts of the impedance,
respectively. The complex-valued impedance is also defined by its _magnitude_
and _phase_ , respectively,

$\displaystyle|Z(\omega)|$
$\displaystyle=\sqrt{Z_{\mathrm{r}}^{2}(\omega)+Z_{\mathrm{j}}^{2}(\omega)}$
(14a) $\displaystyle\angle Z(\omega)$
$\displaystyle=\left\\{\begin{array}[]{ll}\arctan\frac{Z_{\mathrm{j}}(\omega)}{Z_{\mathrm{r}}(\omega)}&\text{for
}Z_{\mathrm{r}}(\omega)\geq 0\\\
\pi+\arctan\frac{Z_{\mathrm{j}}(\omega)}{Z_{\mathrm{r}}(\omega)}&\text{for
}Z_{\mathrm{r}}(\omega)<0.\end{array}\right.$ (14d)

Note that in electrochemistry the phase of the impedance is often denoted as
$\varphi(\omega)$.

The impedance is usually visualised on a Bode plot (Fig. 3 (a-b)) as magnitude
and phase in function of frequency, or on a Nyquist chart (Fig. 3 (c)) as real
versus _negative_ imaginary part, since electrochemical systems are often
capacitative in their electrical behaviour.

Figure 3: Illustration of classical impedance data at different operating
points for a Samsung 48X Li-ion battery (see Section 10). The operating points
in the temperature-SOC plane are shown in (d), while the corresponding
impedances are shown as Bode plot in (a-b) and as Nyquist chart in (c). The
impedance data at the selected frequencies is indicated by dots, which are
connected with straight lines.

In the _potentiostatic_ case, the impedance can still be computed by (12),
since the admittance $Y(\omega)$ is defined as follows,

$\displaystyle Y(\omega)=\frac{I(\omega)}{V(\omega)}=\frac{1}{Z(\omega)}.$
(15)

#### Kramers-Kronig

The conformity of the required constraints of linearity (3) and stationarity
(4) can be validated through the Kramers-Kronig transformation [5, 93, 94,
95]. This transformation states that there is an analytical relation between
the real and imaginary parts of the impedance,

$\displaystyle Z_{\mathrm{r}}(\omega)$
$\displaystyle=Z_{\mathrm{r}}(\infty)+\frac{2}{\pi}\int_{0}^{\infty}\frac{xZ_{\mathrm{j}}(x)-\omega
Z_{\mathrm{j}}(\omega)}{x^{2}-\omega^{2}}\mathrm{d}x$ (16a) $\displaystyle
Z_{\mathrm{j}}(\omega)$
$\displaystyle=-\frac{2\omega}{\pi}\int_{0}^{\infty}\frac{Z_{\mathrm{r}}(x)-Z_{\mathrm{r}}(\omega)}{x^{2}-\omega^{2}}\mathrm{d}x.$
(16b)

In theory, when the imaginary impedance computed from the measured real part
coincides well with the measured imaginary part, or vice versa, the required
constraints are assumed to be satisfied. However, in practice, (16) is
difficult to implement since the integral needs continuous impedance data
(while only discrete data is available) and goes from DC to infinity (while
data is only available in a certain frequency band). Moreover, measurement
noise is not accounted for. Hence, an alternative approach is to fit an
equivalent circuit model based on solely the real, or imaginary, part of the
impedance data. When the model coincides well with the measured complex
impedance data, the Kramer-Kronig relation is assumed to be satisfied [96,
97].

#### Models at operating points

It is very important to stress that the impedance $Z(\omega)$ measured with
classical EIS is dependent on the local operating point (in the sense of a
Taylor expansion) at which the experiments have been performed. The operating
point is defined by the OCV value and the constant values of the external
parameters $p(t)$. As an example, Fig. 3 (d) shows operating points of a Li-
ion battery depending on the state-of-charge (SOC) and temperature. The
measured classical impedances $Z(\omega)$ of a Samsung 48X Li-ion battery (see
Section 10) at these operating points are shown as Bode (a-b) and Nyquist (c)
plots. We observe that the impedance depends on the SOC and temperature, and
that low SOC and low temperature exhibit higher impedance values. Furthermore,
for batteries, a difference in impedance would also be visible for experiments
at the same SOC and temperature, but at different SOH [12, 18].

###  3.2 Excitation

The excitation signal $i(t)$ must be ‘rich’ enough such that the response
$v(t)$ contains the information needed for extracting impedance data. Since
impedance data should be measured at a set of frequencies, it is natural to
use sinusoidal functions as excitations. Historically, EIS was performed with
single-sine excitations. Later, multisine EIS was developed [40]. Both of them
have pros and cons, which are discussed now.

Figure 4: Illustration of the response of an LTI system excited by single-sine
and multisine excitations. The first row has a single-sine excitation at
angular frequency $\omega_{1}$, the second row at $\omega_{2}=2\omega_{1}$,
and the third row has a multisine excitation, which is the sum of the two
single-sines. The OCV is plotted in black.

####  3.2.1 Single-sine excitation

Commonly, single-sine excitations are used for EIS [98, 99]. That is, a
sinusoidal zero-mean current signal with small amplitude $I_{m}$ at a selected
angular frequency $\omega_{m}$ is applied,

$\displaystyle i(t)=I_{m}\cos(\omega_{m}t),$ (17)

and the voltage response (10) is measured,

$\displaystyle
v(t)=\mathrm{OCV}+\underbrace{|Z(\omega_{m})|I_{m}}_{V_{m}}\cos\big{(}\omega_{m}t+\underbrace{\angle
Z(\omega_{m})}_{\psi_{m}}\big{)}.$ (18)

The voltage response is a sinusoidal signal (assuming linearity due to the
small current amplitude) at the same frequency, however, with a different
amplitude $V_{m}$ and phase $\psi_{m}$, superimposed on the OCV. This is
illustrated in Fig. 2 and the two top rows of Fig. 4. The complex-valued
impedance at angular frequency $\omega_{m}$ is computed from the amplitude
scaling and phase shift between current and voltage,

$\displaystyle Z(\omega_{m})=\frac{V_{m}}{I_{m}}e^{j\psi_{m}}.$ (19)

The selected frequencies $\omega_{m}$, $m=1,2,...,M$, are applied
_sequentially_ (i.e. one after the other), usually starting from the highest
frequency and ending at the lowest one. Often the selected frequencies are
logarithmically spaced over multiple decades such as to excite processes
happening at different time-scales. The impedance at each of these frequencies
is computed. Note that it is only possible to apply the sinusoids sequentially
because stationarity (4) is assumed, and, hence, the response is independent
of the time of excitation.

####  3.2.2 Multisine excitation

Instead of applying sinusoidal signals _sequentially_ , it is also possible to
apply the different frequencies _simultaneously_. This is the purpose of a
multisine, where the _sum_ of sinusoidal signals at different frequencies is
applied,

$\displaystyle i(t)=\sum_{m=1}^{M}I_{m}\cos(\omega_{m}t+\phi_{m}).$ (20)

Here, each of the sinusoidal components is given a different phase $\phi_{m}$
such as not to introduce constructive interference (see Section 5). Due to the
linearity constraint (3), the total response yields the sum of each of the
individual responses,

$\displaystyle
v(t)=\mathrm{OCV}+\sum_{m=1}^{M}\underbrace{|Z(\omega_{m})|I_{m}}_{V_{m}}\cos\big{(}\omega_{m}t+\underbrace{\phi_{m}+\angle
Z(\omega_{m})}_{\psi_{m}}\big{)}.$ (21)

The impedance values at the selected angular frequencies $\omega_{m}$ can
still be computed based on the amplitude scaling and the phase shifts of the
sinusoidal components,

$\displaystyle Z(\omega_{m})=\frac{V_{m}}{I_{m}}e^{j(\psi_{m}-\phi_{m})}.$
(22)

An illustration of current and voltage signals for a multisine excitation
under the assumptions of linearity and stationarity is shown in the bottom row
of Fig. 4.

####  3.2.3 The choice of excitation signal

For classical EIS experiments, it is common to use single-sine experiments
since they are known to electrochemists and easy to use with commercially
available potentiostats. However, for broadband experiments, that is,
experiments over a large frequency band, we recommend the use of a multisine.

First of all, since for multisine experiments all frequencies are applied
_simultaneously_ , as opposed to _sequentially_ in single-sine experiments,
the experiment time is shorter [40]. Moreover, in single-sine experiments we
should wait for transients to fade out at each individual frequency, while for
multisine experiments this should only be done once, which also decreases the
experiment time.

The constraints of linearity and stationarity are easily checked for multisine
experiments by looking at the measured current and voltage data in the
frequency domain (see Section 7), while for single-sine experiments one needs
to check the Kramers-Kronig relations. It is noteworthy that for multisine
experiments one must use the frequency domain to detect nonlinearity and
nonstationarity, since it has been empirically shown that multisine
experiments always satisfy the Kramers-Kronig relations [100]. This is studied
in more detail in Section 7. Next, we demonstrate that the most commonly
stated arguments in favor of single-sine experiments over multisine ones can
be debunked.

  1. 1.

The linearity constraint is believed to be easier to impose in the single-sine
setting. The amplitudes such that the response stays linear can be determined
for each of the sines separately. In the multisine case, since a high number
of frequencies are added, the amplitude of each frequency separately should
remain small such that the total multisine signal does not become too high in
magnitude (in a root-mean-square sense) and introduce nonlinearities. However,
finding the optimal excitation amplitudes for multisine experiments can also
be done by detecting nonlinear behaviour. Moreover, the amplitudes can be
dependent on the frequency too.

  2. 2.

The signal-to-noise ratio (SNR) is believed to be better in the single-sine
case. This is since all the power is injected at one frequency, while, in the
multisine case, there is a trade-off between the number of selected
frequencies and the SNR. The more frequencies are excited, the smaller their
amplitudes must be for the response to remain linear, and hence, the smaller
the SNR. However, it is important to take the measurement time into account,
and compare the SNR for single-sine and multisine experiments of the same
duration [26, Section 5.2.2 p. 154]. Since single-sine experiments take a
longer time, we can measure more periods of the multisine during the same
measurement time, resulting in a better SNR.

  3. 3.

For a multisine excitation, the selected frequencies should all be integer
multiples of a fundamental frequency such that the period of the multisine
equals the period of the fundamental sine. These integer multiples are called
_harmonics_. Single-sine experiments do not have this limitation. Accordingly,
multisine experiments have less flexibility in the choice of frequencies at
the low frequency bands. However, this is not a significant issue; one could
perform multiple multisine experiments with slightly different fundamental
frequencies if required.

  4. 4.

The extraction of the impedance in the single-sine case is very intuitive, and
can be handled in the time domain. For multisine excitations, the Fourier
transform is usually used for separating the frequency components, and the
impedance is computed by a ratio of Fourier spectra (see Section 6). Working
in the time domain only, however, is limited because one might not see
leakage, nonlinear distortions and nonstationarity.

Another strong argument in favour of multisine experiments is that if the
system behaves in a nonstationary way, one _should_ use a multisine
excitation, as pointed out in Section 4.2.

In Section 6, we study how classical impedance data can be estimated from
measured current and voltage data.

##  4 Nonlinear and nonstationary impedance models

In this section, impedance models are introduced beyond the very restrictive
constraints of linearity and stationarity. Eliminating the linearity
constraint leads to the Volterra series model, from which the concepts of
nonlinear EIS (NLEIS) and the best linear approximation (BLA) are derived
(Section 4.1). Getting rid of the stationarity constraint leads to the time-
varying impedance model (Section 4.2). Eliminating both the constraints of
linearity and stationarity, we will study the best linear time-varying
approximation (BLTVA) (Section 4.3).

###  4.1 Nonlinear models

The choice of the excitation magnitudes $I_{m}$ in classical EIS measurements
is a compromise between the need to achieve linearity and the need for a
sufficient signal-to-noise ratio [5]. EIS measurements, just as any other
measurements, contain noise. This noise is caused by external disturbances and
the electronics of the measurement device. To improve the quality of impedance
data, one can increase the magnitudes $I_{m}$ in the excitation such that the
voltage response becomes more dominant over the noise level. However, this may
jeopardise the assumption of linearity since the amplitude span over which the
current or voltage is perturbed increases (see Fig. 5). In the case where weak
nonlinear distortions are present in measurements, the estimation of the best
linear approximation (BLA) [36] is a convenient tool. On the other hand, we
could also intentionally introduce nonlinear behaviour in experiments to
investigate additional properties of the electrochemical system, as studied in
NLEIS [30, 31, 33, 34] or intermodulated differential immittance111Impedance
and admittance as a combined concept. spectroscopy [101, 102, 103]. In this
section, we introduce a nonlinear model for EIS measurements through the
Volterra series, derive NLEIS from it, and study the concept of BLA. It is
noteworthy that these are still models at fixed local operating points,
though, but valid over a larger excitation amplitude than classical EIS
measurements.

Figure 5: Illustration of nonlinear time-invariant behaviour in EIS
measurements under a large excitation amplitude. The simulated nonlinearity
comes from the Butler-Volmer equation (7).

####  4.1.1 The Volterra series

Nonlinear time-invariant (NLTI) systems can often be modelled by Volterra
series [27, 28, 104]. Note that this is not the case for all NLTI systems.
Systems with subharmonics or chaotic behaviour, for instance, cannot be
captured by Volterra series. Fortunately, most electrochemical systems can.
Since stationarity is still assumed and we are looking at nonlinear behaviour
at specific local operating points, the current excitation signal should be
zero-mean and the external parameters constant. Mathematically, the general
relation between output voltage and input current (1) may be written as

$\displaystyle v(t)=\mathrm{OCV}+\sum_{n=1}^{n_{\mathrm{max}}}v_{n}(t),$ (23a)
with $n_{\mathrm{max}}$ the order of the nonlinearity, with $\displaystyle
v_{n}(t)=\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}z_{n}(\tau_{1},...,\tau_{n})\prod_{l=1}^{n}i(t-\tau_{l})\mathrm{d}\tau_{l}$
(23b)

being the contribution of the $n$-th order nonlinearity to the voltage output
signal, and $z_{n}(\tau_{1},\ldots,\tau_{n})$ the generalised impulse response
of the $n$-th order nonlinearity.

The Volterra series can be understood as the extension of a Taylor expansion
(8) at a local operating point for dynamical systems. The nonlinear behaviour
is written as a polynomial (instead of linear) operator acting on the
excitation current.

The generalised impedances $Z_{n}$ are defined as the $n$-dimensional Fourier
transform of the generalised impulse responses $z_{n}$:

$\displaystyle
Z_{n}(\omega_{1},...,\omega_{n})=\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}z_{n}(\tau_{1},\ldots,\tau_{n})\prod_{l=1}^{n}e^{-j\omega_{l}\tau_{l}}\mathrm{d}\tau_{l}.$
(24)

For nonlinearity order $n_{\mathrm{max}}=1$, the current-voltage relations of
linear systems (9) and (10) are retrieved from the Volterra series,

$\displaystyle v_{1}(t)$
$\displaystyle=\int_{-\infty}^{\infty}z_{1}(\tau)i(t-\tau)\mathrm{d}\tau=\mathcal{F}^{-1}\\{Z_{1}(\omega)I(\omega)\\}.$
(25)

##### Response to a zero-mean single-sine

Let us now look at the response of an NLTI system, described by a Volterra
series, to a zero-mean single-sine excitation $i(t)=I\cos(\omega t)$. The
contribution of the $n$-th order nonlinearity to the voltage signal yields
[104]

$\displaystyle v_{n}(t)=\sum_{h=0}^{n}|V_{n,h}|\cos\big{(}h\omega t+\angle
V_{n,h}\big{)},$ (26a) with $\displaystyle V_{n,h}$
$\displaystyle=Z_{n,h}(\omega)I^{n}\qquad h>0$ (26b) $\displaystyle
Z_{n,h}(\omega)$
$\displaystyle=\frac{1}{2^{n-1}}\sum_{\\{s_{1},\ldots,s_{n}\\}\in\mathbb{S}_{n,h}}Z_{n}(s_{1}\omega,\ldots,s_{n}\omega).$
(26c)

For harmonic $h=0$, the premultiplying factor should be $1/2^{n}$ instead of
$1/2^{n-1}$. The set $\mathbb{S}_{n,h}$ contains all possible lists
$\\{s_{1},\ldots,s_{n}\\}$, with elements $s_{1,\dots,n}\in\\{-1,1\\}$, such
that $s_{1}+\ldots+s_{n}=h$. When the set is empty, $Z_{n,h}(\omega)=0$. The
values $Z_{n,h}(\omega)$ are called the _nonlinear impedance coefficients_ ,
and come directly from the generalised impedances. A sum of $n$ elements with
values that are either $1$ or $-1$ can never be larger than $n$ and the sum of
an even number of these elements can never be odd, and vice versa. This
translates into,

$\displaystyle Z_{n,h}(\omega)=0$ $\displaystyle\text{for }h>n$ (27a)
$\displaystyle Z_{2n,2h+1}(\omega)=0$ $\displaystyle\forall n,h\in\mathbb{N}$
(27b) $\displaystyle Z_{2n+1,2h}(\omega)=0$ $\displaystyle\forall
n,h\in\mathbb{N}.$ (27c)

These seemingly complicated mathematics are better understood by looking at a
few special cases. For the linear term ($n=1$) in the Volterra series we find
that $V_{1,0}=0$ and $\mathbb{S}_{1,1}=\\{1\\}$, such that
$V_{1,1}=Z_{1}(\omega)I$. Hence, the voltage response yields

$\displaystyle v_{1}(t)$
$\displaystyle=\underbrace{|Z_{1}(\omega)|I}_{|V_{1,1}|}\cos\big{(}\omega
t+\underbrace{\angle Z_{1}(\omega)}_{\angle V_{1,1}}\big{)}.$ (28)

The response is present at the same frequency as the excitation, which is in
accordance with the expected linear output (18). For nonlinear systems, that
is $n_{\mathrm{max}}\geq 2$, we notice from (26a) that spectral content may be
present at _integer_ multiples of the excited frequency. Considering (23) and
(27) with purely quadratic ($n=2$) and purely cubic ($n=3$) nonlinearities, we
find, respectively,

$\displaystyle v_{2}(t)$ $\displaystyle=|V_{2,0}|+|V_{2,2}|\cos\big{(}2\omega
t+\angle V_{2,2}\big{)}$ $\displaystyle v_{3}(t)$
$\displaystyle=|V_{3,1}|\cos\big{(}\omega t+\angle
V_{3,1}\big{)}+|V_{3,3}|\cos\big{(}3\omega t+\angle V_{3,3}\big{)}.$ (29)

The expressions for the different $V_{x,y}$ terms are given in B. We remark
that quadratic nonlinearities introduce spectral content at the _even_ integer
multiples of the excited frequency smaller or equal to two, and that cubic
nonlinearities introduce spectral content at the _odd_ integer multiples of
the excited frequency smaller or equal to three. These results can be
generalised to higher order even and odd nonlinearities,

$\displaystyle v_{2n}(t)$
$\displaystyle=\sum_{h=0}^{n}|V_{2n,2h}|\cos\big{(}2h\omega t+\angle
V_{2n,2h}\big{)}$ $\displaystyle v_{2n+1}(t)$
$\displaystyle=\sum_{h=0}^{n}|V_{2n+1,2h+1}|\cos\big{(}(2h+1)\omega t+\angle
V_{2n+1,2h+1}\big{)}.$ (30)

Accordingly, the even and odd degree terms in the Volterra series introduce
frequency content at, respectively, even and odd integer multiples of the
excited frequency. Even nonlinear functions can be captured by the sum of even
degree monomials of the Volterra series whereas odd nonlinear functions can be
captured by odd degree monomials. Accordingly, even nonlinear behaviour is
present at even integer multiples of the excited frequency, and odd nonlinear
behaviour at odd multiples.

In contrast with the small excitation magnitude applied to the Butler-Volmer
equation (Fig. 2), where the response is only present at the excited frequency
$f_{1}$, for a larger excitation magnitude (Fig. 5), spectral content may also
be present at integer multiples of the excited frequency. For this particular
illustrative example, only odd nonlinear distortions are present, since we
fixed $\alpha_{a}=\alpha_{c}=0.5$ and therefore the Butler-Volmer equation (6)
is an odd function around $0$. However, in practice, when
$\alpha_{a}\neq\alpha_{c}$ the Butler-Volmer equation is neither even nor odd,
and hence, both even and odd nonlinear distortions are present.

Writing out the analytical expression of the response of an NLTI system
described by a Volterra series to a _multisine_ excitation is a more
complicated matter. The mathematical details are omitted for this review
paper, however, they are given in Lang and Billings [104]. Fortunately, when
the multisine consists of excited frequencies at integer multiples of a
fundamental frequency, the observations made above are still valid.
Nonlinearities will still be present at integer multiples of this fundamental
frequency. However, the distinction between even and odd nonlinear distortions
can only be made when only odd integer multiples of the fundamental frequency
are excited in the multisine. This is discussed in Section 7.2.

####  4.1.2 Nonlinear EIS

It can be shown that the total response of the Volterra series of infinite
order ($n_{\mathrm{max}}\rightarrow\infty$), excited by a sinusoidal signal at
frequency $\omega$, $i(t)=I\cos(\omega t)$, introduces spectral content at all
the integer multiples of the excited frequency $\omega$, that is,

$\displaystyle v(t)=\mathrm{OCV}+\sum_{h=0}^{\infty}|V_{h}|\cos(h\omega
t+\angle V_{h}),$ (31a) with $\displaystyle V_{h}$
$\displaystyle=\sum_{n=1}^{\infty}V_{n,h}=\sum_{n=h}^{\infty}Z_{n,h}(\omega)I^{n}$
(31b) $\displaystyle=\sum_{r=0}^{\infty}Z_{h+2r,h}(\omega)I^{h+2r}.$ (31c)

In these equations, (23), (26) and (27) have been exploited. Only even order
nonlinear impedance coefficients larger or equal to $h$ introduce spectral
content at an even harmonic $h$, and vice versa for odd harmonics.

Nonlinear EIS [105, 31, 34] aims at measuring the _leading order_ nonlinear
impedance coefficients $Z_{h,h}(\omega)$. These are defined from (31) as [34],

$\displaystyle Z_{h,h}(\omega)=\lim_{I\rightarrow 0}\frac{V_{h}}{I^{h}}.$ (32)

Note that the unit of the $h$-th leading order nonlinear coefficient
$Z_{h,h}(\omega)$ is $\Omega/\text{A}^{h-1}$. In practice, the leading order
nonlinear coefficients are measured as,

$\displaystyle\hat{Z}_{h,h}(\omega)=\frac{V_{h}}{I^{h}}=Z_{h,h}(\omega)+\underbrace{\sum_{r=1}^{\infty}Z_{h+2r,h}(\omega)I^{2r}}_{\text{choose
$I$ such that negligible}},$ (33)

where the excitation amplitude $I$ is chosen large enough such that the $h$-th
harmonic $V_{h}$ is visible, but also small enough such that the higher order
contributions $Z_{h+2r,h}(\omega)I^{r}$, $r\geq 1$, are negligible. Note that
under these measurement conditions, $Z_{1,1}(\omega)$ is the regular impedance
$Z(\omega)$. The second leading order nonlinear impedance coefficient was
measured for Li-ion batteries in [31] and [34]. Recall that the leading order
nonlinear impedance coefficients are still models at a specific local
operating point.

####  4.1.3 The best linear approximation

For reasons of simplicity, we may prefer to work with linear models. In
Schoukens et al. [36], NLTI systems are modelled as so-called best linear
approximations, plus a ‘nonlinear noise source’ generating nonlinear
distortions $v_{\mathrm{s}}(t)$. Hence, in this context

$\displaystyle
v(t)=\mathrm{OCV}+\underbrace{\mathcal{F}^{-1}\\{Z(\omega)I(\omega)\\}}_{v_{\mathrm{BLA}}(t)}+v_{\mathrm{s}}(t),$
(34)

where $v_{\mathrm{BLA}}(t)$ stands for the response of the BLA $Z(\omega)$.
The BLA is the ‘best’ linear model in the sense that it minimises the
nonlinear distortions in least square sense,

$\displaystyle
Z=\arg\min_{Z^{\prime}}\mathbb{E}\Biggl{\\{}\int_{-\infty}^{\infty}|v_{\mathrm{s}}(t)|^{2}\mathrm{d}t\Biggr{\\}},$
(35)

with the expected value $\mathbb{E}\\{\cdot\\}$ taken over different
realisation if the excitation signal [106]. Using Parceval’s theorem, the BLA
can also be defined in the frequency domain,

$\displaystyle Z(\omega)$
$\displaystyle=\arg\min_{Z^{\prime}(\omega)}\mathbb{E}\Bigl{\\{}|V(\omega)-Z^{\prime}(\omega)I(\omega)|^{2}\Bigr{\\}}.$
(36)

It follows that the nonlinear distortions are uncorrelated with, but not
independent of, the excitation signal. Hence, for a zero-mean single-sine
excitation, the response of the BLA is located at the excited frequency, while
the nonlinear distortions are present at the remaining integer multiples of
the excited frequency,

$\displaystyle v_{\mathrm{BLA}}(t)$ $\displaystyle=|V_{1}|\cos(\omega t+\angle
V_{1})$ (37a) $\displaystyle v_{\mathrm{s}}(t)$
$\displaystyle=|V_{0}|+\sum_{h=2}^{\infty}|V_{h}|\cos(h\omega t+\angle
V_{h}),$ (37b)

with the coefficients $V_{h}$ defined in (31). An illustration of this
decomposition into a linear response and purely nonlinear response is depicted
in Fig. 6. The linear response is the one at the excited frequency, while the
nonlinear response consists of the remaining frequencies. This BLA method has
the major advantage that a linear model can still be justified when the
nonlinear distortions are sufficiently small. Accordingly, it still makes
sense to measure an impedance from data which shows nonlinear distortions, as
long as these nonlinear distortions are small enough for the intended
application. Therefore, it is important to detect and quantify the nonlinear
distortions in measurements.

Figure 6: Decomposition of the nonlinear time-invariant response of Fig. 5
into a linear response and purely nonlinear response. Left graph: frequency
domain, right graph: time domain.

The BLA under a single-sine excitation is defined from (31) as,

$\displaystyle
Z(\omega)=\frac{V_{1}}{I}=\sum_{r=0}^{\infty}Z_{1+2r,1}(\omega)I^{2r},$ (38)

and hence, _the BLA depends on the amplitude of the excitation_. This is
illustrated for the static case in Fig. 7, the linearisation depends on the
span over which the Butler-Volmer equation is linearised. The BLA here is the
slope between voltage and current, which is clearly different for the black
and grey lines. The reasoning in the dynamic case is similar.

Note that for nonlinear systems the first generalised impedance and the first
leading order nonlinear impedance coefficients are equal, but the BLA
impedance is different: $Z_{1}(\omega)=Z_{1,1}(\omega)\neq Z(\omega)$. This is
because the higher order odd polynomials in the volterra series also contain a
linear part, leading to a spectral contribution at the excited frequency that
contributes to the BLA $Z(\omega)$.

Figure 7: The linearisation of a static nonlinear function depends on the span
it is linearised over. The red line represents the Butler-Volmer equation (7).
Linearising it over the black span gives a different BLA (the slope) than
linearising it around the grey span.

###  4.2 Nonstationary models

The stationarity constraint in classical EIS experiments is very restrictive.
The operating point, and hence also the external parameters, should be
constant during an experiment. (NL)EIS can only be performed on systems in
steady-state, resulting in models at fixed operating points. Accordingly, if
we want impedance data over various operating conditions, for instance at
different OCV values as in Fig. 3, we have to separately drive the system to
each of these operating conditions and wait for steady-state, which is very
time-consuming. Moreover, sometimes the system never reaches steady-state due
to inherently nonstationary behaviour caused by changes in thermodynamic
states and kinetically slow side processes (such as the self discharge of
energy storage systems). Furthermore, it is of great interest to study
electrochemical systems during _operation_. Examples include the formation of
film layers during anodising, the electrorefining of copper, and the charging,
discharging, and relaxation of batteries. For these examples, classical EIS or
NLEIS can only be performed as a perturbation around a rest condition, while
the evolution of the impedance during operation contains important
information. Unfortunately, this information cannot be gathered by classical
or stationary nonlinear EIS.

Nonstationarity can occur for two reasons, and they can happen simultaneously.
The first cause is that the external parameters $p(t)$ vary during the
experiment. The other cause is that the system is excited in such a way that
it does not remain in steady-state during the experiment. This happens when
superimposing a conventional excitation $i_{\text{exc}}(t)$ (see Section 3.2)
on a slow signal $i_{0}(t)$, driving the system in operating conditions,

$\displaystyle i(t)=i_{0}(t)+i_{\text{exc}}(t).$ (39)

For a battery, time-variation, for instance, occurs when the excitation is a
multisine superimposed on a constant offset $i_{0}$ that (dis)charges the
battery [59], or when a zero-mean excitation is applied right after charging
or discharging to study the relaxation behaviour [61].

####  4.2.1 The time-varying impedance

Considering the two sources of nonstationarity, the linear voltage response
can be modelled as,

$\displaystyle
v(t)=v_{0}(t)+\int_{-\infty}^{\infty}z(\tau,t)i_{\mathrm{exc}}(\tau)\mathrm{d}\tau.$
(40)

Here $v_{0}(t)$ represents a drift signal. In the battery example, this would
be the voltage slowly increasing as the battery is charging due to a positive
constant current. The time-variation due to the external parameters $p(t)$
and/or excitation trajectory $i_{0}(t)$ are simultaneously captured by a two-
dimensional impulse response. This time-varying impulse response $z(\tau,t)$
was introduced by Zadeh [37] in 1950 for modelling LTV systems. It is a
natural extension of (9), where now the impulse response function explicitly
depends on the excitation time.

It is also shown in Battistel et al. [107] that nonstationarity (40) appears
from an NLTI system when the response is linearised along a time-varying
trajectory. This is detailed in C.

Similarly to (10), the time-varying impedance $Z(\omega,t)$ appears when
transforming (40) into the frequency domain [37],

$\displaystyle
v(t)=v_{0}(t)+\mathcal{F}^{-1}\\{Z(\omega,t)I_{\mathrm{exc}}(\omega)\\},$ (41)

with the time-varying impedance defined as,

$\displaystyle
Z(\omega,t)=\int_{-\infty}^{\infty}z(t-\tau,t)e^{-j\omega\tau}\mathrm{d}\tau.$
(42)

Accordingly, when a single-sine excitation (17) is applied, superimposed on a
slowly varying trajectory $i_{0}(t)$, the voltage response yields,

$\displaystyle
v(t)=v_{0}(t)+|Z(\omega_{m},t)|I_{m}\cos\left(\omega_{m}t+\phi_{m}+\angle
Z\left(\omega_{m},t\right)\right).$ (43)

The current excitation is, hence, modulated in amplitude and phase by
$Z(\omega_{m},t)$. Since linearity is still assumed, the response to a
multisine excitation (20) is simply the sum of the responses to each separate
sinusoidal signal,

$v(t)=v_{0}(t)\>+\>\\\
\sum_{m=1}^{M}|Z(\omega_{m},t)|I_{m}\cos\left(\omega_{m}t+\phi_{m}+\angle
Z\left(\omega_{m},t\right)\right).$ (44)

Figure 8: Illustration of the response of a Li-ion battery (which can be
modelled by a Volterra series) to zero-mean and nonzero-mean excitations with
small amplitudes. The slow parts $i_{0}(t)$ and $v_{0}(t)$ are shown in black.
The positive mean value of the current charges the battery, and, hence, the
voltage increases.

An illustration of the response of a Li-ion battery under zero-mean and
nonzero-mean small amplitude excitation is shown in Fig. 8. For the zero-mean
excitation, the response can be modelled as the OCV plus an LTI response. For
the nonzero-mean excitation, nonstationarity is introduced due to the
constant-current charging, which also might cause external parameters such as
the temperature to change. Accordingly, the response can be modelled as a
drift signal plus the response of an LTV system.

####  4.2.2 Models along an operating trajectory

By choosing the slow trajectory $i_{0}(t)$ and/or varying external parameters
during the experiment, we obtain a linear model for the impedance along an
operating _trajectory_ instead of an operating _point_. This trajectory is the
drift of the system due to external effects during the measurement. One can,
hence, obtain more global models than with classical EIS. This is illustrated
in Fig. 9, where a battery with capacity $C$ =
$4.8\text{\,}\mathrm{A}\mathrm{h}$ is charged using a $C/2$ current, that is,
$i_{0}(t)$ = $2.4\text{\,}\mathrm{A}$. Note that here the external temperature
(a parameter), measured at the battery surface, changes due to thermal
dynamics. The SOC also changes, however, this is not an external parameter,
since it depends on the current excitation.

Figure 9: Illustration of time-varying impedance data of a Samsung 48X battery
along an operating trajectory. The operating trajectory is caused by a
charging current $i_{0}(t)=2.4$ A applied to a $4.8\,\mathrm{Ah}$ battery
placed in a thermal chamber at $5^{\circ}$C. The time-varying impedance
$Z(\omega,t)$ along the trajectory in the temperature-SOC plane (d) is shown
as Bode plot (a-b) and Nyquist plot (c).

####  4.2.3 The importance of multisine excitation

A multisine excitation is mandatory for accurate estimation of time-varying
impedance. Since the system is changing over time during the experiment, it
would be illogical to apply single-sines _sequentially_ , since then the
impedance at each selected frequency would be computed for a different section
of the operating trajectory. The advantage of using a multisine excitation is
that many frequencies are excited _simultaneously_ , and we can obtain the
impedance at all the selected frequencies over the _entire_ operating
trajectory. In Section 9, we study how time-varying impedance data can be
estimated from measured current and voltage data under multisine excitation.

###  4.3 Nonlinear and nonstationary models

When a system is excited whilst not in steady-state, and with large excitation
amplitudes $I_{m}$, nonstationary and nonlinear behaviour may happen
_simultaneously_. The system is then denoted as nonlinear time-varying (NLTV).
In this case, the time-varying Volterra series with a superimposed drift
signal provides a general model for the response,

$\displaystyle
v(t)=v_{0}(t)+\sum_{n=1}^{n_{\mathrm{max}}}\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}z_{n}(\tau_{1},...,\tau_{n},t)\prod_{l=1}^{n}i_{\mathrm{exc}}(\tau_{l})\mathrm{d}\tau_{l}.$
(45)

Ideally, from such a time-varying Volterra series, we could measure time-
varying leading order nonlinear impedance functions $Z_{h,h}(\omega,t)$. This
would provide models over a large excitation amplitude and along a time-
varying trajectory. However, this has, to the best of our knowledge, not been
studied yet.

In Hallemans et al. [49] we have studied the extension of the concept of BLA
to BLTVA, that is, the best linear time-varying approximation. In this
framework, the relation between current and voltage of an NLTV system is
modelled as,

$\displaystyle
v(t)=v_{0}(t)+\underbrace{\mathcal{F}^{-1}\\{Z(\omega,t)I_{\mathrm{exc}}(\omega)\\}}_{v_{\mathrm{BLTVA}}(t)}+v_{\mathrm{s}}(t),$
(46)

with $Z(\omega,t)$ the BLTVA and $v_{\mathrm{s}}(t)$ the time-varying
nonlinear distortions. The BLTVA is a promising tool to monitor
electrochemical system impedance during operation [63, 59] (see Section 9.3).

##  5 Measuring current and voltage data

In practical settings, we cannot measure continuous-time signals $i(t)$ and
$v(t)$ over infinite periods. Instead, sampled current and voltage data over
finite periods should be collected, denoted $i(n)$ and $v(n)$ where $n$ is the
sample number. To obtain these, we apply an excitation through a potentiostat
and measure sampled and windowed current and voltage data. For extracting
time-varying impedance from this data, a multisine excitation is recommended,
as used for odd random phase (ORP) EIS [59] and dynamic multi-frequency
analysis [107].

##### Design of excitation signal

Different kinds of excitation signals are ‘rich’ enough to estimate classical
impedance; among others, single-sines, multisines, and white noise. For
obtaining stationary nonlinear impedance estimates, a single-sine excitation
should be used. For obtaining time-varying impedance data, a multisine should
be used. Since the single-sine excitation is a special case of the multisine,
we focus on the latter.

A multisine with period $T_{p}$ superimposed on a time-varying trajectory
$i_{0}(t)$ is given by

$\displaystyle
i(t)=i_{0}(t)+\underbrace{\sum_{m=1}^{M}I_{m}\cos\Big{(}\frac{2\pi
h_{m}}{T_{p}}t+\phi_{m}\Big{)}}_{i_{\mathrm{exc}}(t)}.$ (47)

The trajectory $i_{0}(t)$ is user defined—for example this could be a charging
or discharging current for a battery, or a chronoamperometry trajectory. Note
that to obtain classical or stationary NLEIS data, this trajectory should be
zero. The set of excited harmonics is defined as
$\mathbb{H}_{\text{exc}}=\\{h_{1},h_{2},...,h_{M}\\}$, and the excited angular
frequencies are accordingly $\omega_{h_{m}}=2\pi h_{m}/T_{p}$, with $T_{p}$
the period of the multisine. The harmonic numbers should be integers,
$h_{m}\in\mathbb{N}$, such that all sinusoidal signals fit an integer number
of times in the period $T_{p}$. Note that the lowest frequency in the
multisine ($f_{1}=1/T_{p}$) is inversely proportional to the period of the
multisine. For a single-sine excitation, only $f_{1}$ is excited. In our
definition, the natural numbers $\mathbb{N}$ do not include zero, while the
set $\mathbb{N}_{0}$ does include zero. The amplitudes are selected by the
user, depending on the application. The phases are most of the time chosen
such as to minimise the crest factor of the overall multisine, that is, as not
to introduce constructive interference when adding sines to each other.
Different approaches can be used for this, including random phases picked from
a uniform distribution in $[0,2\pi)$, and mathematical optimisation techniques
such as the Schröder phase or DFT-based iterative algorithms to minimise the
crest factor [108, 109, 110, 111, 112, 57].

The name _odd random phase_ electrochemical impedance spectroscopy denotes
impedance measurements under multisine excitation with only odd harmonics
excited ($h_{m}\in 2\mathbb{N}_{0}+1$), with random phases. This excitation
signal was introduced by Hubin and Pintelon et al. [40, 41]. In Section 7, we
show that inherent nonlinearity and nonstationarity in electrochemical systems
can easily be detected using an ORP multisine excitation, and we discuss why
it is advantageous to only excite _odd_ harmonics.

##### Windowing

We cannot measure signals for an infinitely long time, but only for a certain
period $t\in[0,T)$. The measurement time $T$ is chosen to measure a certain
ongoing reaction, for instance a charge cycle of a Li-ion battery. To avoid
spectral leakage in the frequency domain, an _integer_ number of periods of a
multisine excitation should be measured, that is, $T=PT_{p}$ with
$P\in\mathbb{N}$ [26, Section 2.2.3, p 40]222Note that measuring an integer
number of periods is only a requirement when the impedance estimation is
performed in the frequency domain, which is the case for multisine
experiments, but not necessarily for single-sine experiments.. This is not
always possible, but it is strongly recommended. Moreover, for obtaining NLEIS
or time-varying impedance data, measuring an integer number of periods is a
requirement.

##### Sampling

Only a sampled representation of the continuous signal can be recorded, at a
sampling frequency $f_{s}$. Following the Shannon-Nyquist sampling theorem,
this sampling frequency should be greater than twice the highest frequency in
the measurements to avoid spectral aliasing. The sampling period is
$T_{s}=1/f_{s}$. It is important that the data is uniformly sampled.

##### Measuring the data

The sampled and windowed multisine current data is applied to an
electrochemical device using a potentiostat. The potentiostat uses a digital-
to-analog converter (DAC) to transform the generated time-series to a
continuous signal. User-defined excitation is not always available in
commercial potentiostats, but user-defined excitation is essential for the
techniques in this article. Accordingly, multiple periods of the multisine
excitation signal can be applied, and the potentiostat then measures the
actual current and the voltage, which are also windowed and sampled.

The collected data can be written as follows,

$\displaystyle\mathcal{D}_{\text{time}}=\left\\{\begin{array}[]{ll}&[i(0),i(1),\ldots,i(N-1)]\\\
&[v(0),v(1),\ldots,v(N-1)]\end{array}\right\\},$ (50)

where $x(n)$ is shorthand notation for $x(nT_{s})$, $x=i,v$. The number of
samples is $N=Tf_{s}$.

##### frequency domain data

Within the constraints of LTI systems, the impedance is defined as the ratio
of the Fourier transforms of voltage and current (12). Hence, it would be
appropriate to directly compute the impedance in the frequency domain.
However, the Fourier transform acts on continuous signals, while only
discrete-time measurements (50) can be collected from potentiostats.
Fortunately, the spectrum of time-series can be computed by replacing the
Fourier integral by a discrete sum. This is called the discrete Fourier
transform333In 1965, Cooley and Tukey designed a highly-efficient algorithm to
compute the DFT, which rapidly popularised frequency domain signal processing.
This algorithm became known as the ‘fast Fourier transform’ (FFT) [113] and is
still used to date. (DFT),

$\displaystyle X(k)=\frac{1}{N}\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}\quad
k=0,1,...,N-1.$ (51)

Here, $x(n)$ is again shorthand for $x(nT_{s})$ and $x=i,v$. The DFT index $k$
corresponds to the angular frequency $\omega_{k}=2\pi k/T$ or frequency
$f_{k}=k/T$. The frequency domain current and voltage data yields,

$\displaystyle\mathcal{D}_{\text{freq}}=\left\\{\begin{array}[]{ll}&[I(0),I(1),\ldots,I(N-1)]\\\
&[V(0),V(1),\ldots,V(N-1)]\end{array}\right\\}.$ (54)

When periodic time domain data is measured for an integer number of periods
and sampled satisfying Shannon-Nyquist’s theorem, the DFT coincides with the
Fourier transform, evaluated on the DFT grid $f_{k}=k/T$ with $k=0,1,...,N-1$
[26, Section 2.2, p 34].

An illustration of one period of a windowed and sampled odd random phase
multisine signal in the time and frequency domain, with the $7$-th harmonic
unexcited for detection of odd nonlinear distortions, is shown in the left
plots of Fig. 10. For the time domain plots, the full line is the continuous
signal, while the dots are the sampled data. For the frequency domain plots,
the DFT grid is indicated by vertical red bars on the frequency axis. Since
DFT lines are intentionally left open at even harmonics and the odd harmonic
$7$, it is possible to detect, quantify and classify nonlinear distortions in
the response to one period of the excitation. However, if the system is also
nonstationary, it is not possible to distinguish between nonlinear and
nonstationary behaviour when measuring only one period.

Figure 10: A sampled and windowed odd random-phase multisine in time and
frequency domain. Top row: continuous signal (full line) and sampled data
(blue dots). Bottom row: DFT of the sampled data, with the DFT grid indicated
by red (original grid) and black (grid for $P=3$ measured periods) vertical
bars on the frequency axis.

##### Increasing the frequency resolution

Measuring more than one period increases the _frequency resolution_ , that is,
the distance between two DFT lines $f_{s}/N=1/T$ becomes smaller. This is
illustrated in the right plots of Fig. 10. For one measured period, the DFT
grid is indicated in red, for three measured periods, in black. Nonlinearities
in the spectrum of the response can _only_ be present at the DFT lines
indicated in red, since these are the integer multiples of the fundamental
excitation frequency $1/T_{p}$. The effect of nonstationarities can be present
at all DFT lines, however. Accordingly, measuring a large number of periods is
commonly used to distinguish between nonlinear and nonstationary behaviour,
and the nonstationary behaviour can be modelled from the response spectrum.

##### Influence of measurement noise

We usually assume that noise in the measured data, $i_{\text{meas}}(n)$ and
$v_{\text{meas}}(n)$, is additive,

$\displaystyle x_{\text{meas}}(n)$ $\displaystyle=x(n)+\mathrm{n}_{x}(n)\qquad
x=i,v,$ (55)

and the noise time-series $\mathrm{n}_{x}(n)$ is iid (independent and
identically distributed), zero-mean and Gaussian:

$\displaystyle\mathrm{n}_{x}(n)\sim\mathcal{N}(0,\sigma^{2}_{\mathrm{n}_{x}})\qquad
x=i,v.$ (56)

Here, $\sigma^{2}_{\mathrm{n}_{i}}$ and $\sigma^{2}_{\mathrm{n}_{v}}$ are the
noise variances on the current and voltage, respectively. The DFTs of the
noise time series, $\mathrm{N}_{I}(k)$ and $\mathrm{N}_{V}(k)$ are circular
complex Gaussian444A two-dimensional Gaussian distribution on the complex
plane. distributed. However, the variances scale inversely with the number of
samples $N$,

$\displaystyle\mathrm{N}_{X}(k)\sim\mathcal{N}_{c}\Big{(}0,\underbrace{\frac{\sigma^{2}_{\mathrm{n}_{x}}}{N}}_{\sigma^{2}_{\mathrm{N}_{X}}(k)}\Big{)}\qquad
x=i,v.$ (57)

We can interpret this result as follows: by measuring a higher number of
samples $N$, the number of DFT lines increases, and, hence, the noise is
distributed over more DFT lines, such that the variance of the circular
complex distributed noise at each DFT line decreases. As a consequence, we can
show that the frequency domain SNR increases with the square root of the
number of measured samples,

$\displaystyle\mathrm{SNR}_{X}(k)=\sqrt{\frac{|X(k)|^{2}}{\sigma^{2}_{\mathrm{N}_{X}}(k)}}=\sqrt{N}\frac{|X(k)|}{\sigma_{n_{x}}}.$
(58)

Accordingly, by measuring a higher number of samples $N$, we increase the SNR.
In practice, the SNR is improved by increasing the number of measured periods
$P$ of the data.

##  6 Classical frequency domain impedance estimation

In the early seventies, Creason and Smith [114, 115, 116] adopted the recently
discovered FFT for classical EIS directly from frequency domain. These
techniques, referred to as FFT-EIS, measure the impedance starting from DFT
data (54). An important question is the choice of the excitation signal. As
explained earlier, a zero-mean excitation is needed for stationarity. In
[116], different zero-mean excitation signals are studied: periodic
excitations, transient inputs, and band-limited white noise. The conclusion is
drawn that periodic excitation, with excited frequencies lying on the DFT
grid, are superior. However, it is not always possible to apply periodic
excitation and measure over an integer number of periods. Think for instance
about a low-cost measurement apparatus for a BMS where only short and fixed-
length data records can be handled. Hence, estimating techniques are also
needed for random excitations. In 1975, Blanc et al. [117] proposed the so-
called ‘pseudo-white noise’ technique (or ‘méthode du bruit blanc’, since this
article was written in French). Pseudo-white noise may also be used as an
excitation, and the impedance is computed as the ratio of the cross- and auto
power spectra in the frequency domain. Howey et al. [88] also used the ratio
of the cross- and auto power spectra in the frequency domain to perform fast
impedance measurements on a self-made low-cost excitation and measurement
system for batteries that could be used in a BMS. The same issue of the choice
of excitation signal for frequency domain system identification is studied in
[108] and in Pintelon and Schoukens’ book [26, Chapter 2]. We now
mathematically formalise the impedance estimation for periodic and random
excitation signals.

##### Periodic excitation

When measuring an integer number of periods $P$ in steady-state under a
periodic excitation, for instance a multisine, the DFT is exactly a sampled
version of the continuous Fourier transform [26, Section 2.2, p 34]. Hence,
the impedance can simply be computed from (12) as the ratio of the voltage and
current spectra,

$\displaystyle\hat{Z}(\omega_{k})=\frac{V(k)}{I(k)}$
$\displaystyle\omega_{k}=\frac{2\pi k}{T},\ k\in P\mathbb{H}_{\text{exc}},$
(59)

with $\mathbb{H}_{\text{exc}}$ the excited harmonics of the excitation signal.
A measurement example is shown in Fig. 11. Since the OCV in the voltage is
only present at DC (zero frequency), it has no influence on the estimate of
the impedance at positive frequencies. Note that in perfect LTI conditions,
only measurement noise causes an uncertainty on the estimate (59). Increasing
the number of periods $P$ generates an averaging effect, which reduces this
uncertainty [26, Section 2.4, p 44]. Periodic excitation is recommended when
possible.

Figure 11: Example of classical impedance estimation for a Li-ion battery at
$10$% SOC and $25^{\circ}$C under a periodic multisine excitation measured for
$P=10$ periods. The excited lines, $k\in P\mathbb{H}_{\text{exc}}$, are
indicated with vertical lines with large dots and crosses. The remaining DFT
lines (small dots) contain noise and a small drift signal.

##### Random excitations

For random excitations, such as pseudo-white noise, pseudo-random binary
sequences (PRBS) or multisine excitations not measured for an integer number
of periods, the DFT does not correspond to the continuous Fourier transform
anymore due to transients. These transients can be reduced by using windows,
e.g. the Hann window. For random excitation, it is recommended to measure the
impedance by the ratio of the cross- $\hat{S}_{VI}(k)$ and auto-spectra
$\hat{S}_{II}(k)$ [26, Section 2.6, p. 54],

$\displaystyle\hat{Z}(\omega_{k})=\frac{\hat{S}_{VI}(k)}{\hat{S}_{II}(k)}=\frac{\sum_{m=1}^{M}V_{[m]}(k)I^{*}_{[m]}(k)}{\sum_{m=1}^{M}I_{[m]}(k)I^{*}_{[m]}(k)},$
(60)

where the superscript ∗ stands for the complex conjugate, and the subscript
[m] stands for different experiments. Averaging over many different
experiments reduces the transient and avoids divisions by zero. Note that this
only makes sense in stationary conditions where multiple experiments can be
gathered.

As a special case, Blanc [117] chooses the pseudo-white noise excitation
$i(t)$ such that $I(k)I^{*}(k)=1$ in the frequency band of interest.
Accordingly, only the numerator of (60) needs to be computed to estimate the
impedance. Note that these cross- and auto spectra correspond to correlations
in the time domain [118], as also studied by Blanc [117].

Other, more involved, frequency domain techniques are possible for measuring
the impedance using random excitation [119, 120]. Here, local parametric
modelling and Gaussian process regression are used to separate the impedance
from the transients. Such techniques are promising for low-cost experiments
where periodic excitation signals cannot be applied. However, explaining these
in detail goes beyond the scope of this article.

##  7 Detection of nonlinearity and nonstationarity

It is empirically shown in You et al. [100] that classical impedance data
$\hat{Z}(\omega_{k})$ obtained from broadband excitation using frequency
domain techniques (as studied in Section 6) always satisfies the Kramers-
Kronig relations. Hence, the latter cannot assess to which extent measured
multisine data satisfies the assumptions of linearity and stationarity, and
another tool is required for this purpose. Inherently nonlinear and
nonstationary behaviour of electrochemical systems is easily detected by
applying a zero-mean ORP excitation (47) (with $i_{0}(t)=0$) and studying the
measured frequency domain data (54) [121, 40, 42, 49]. We recommend the use of
this tool, which we now illustrate on a simplistic example that demonstrates
its advantages.

###  7.1 An example with ORP multisine excitation

Consider an odd random phase multisine with excited harmonics
$\mathbb{H}_{\mathrm{exc}}=\\{1,3,7\\}$, that is, excited frequencies
$1/T_{p}$, $3/T_{p}$ and $7/T_{p}$. A measurement is performed for a duration
$T=10T_{p}$, i.e., $P=10$ periods are measured. Based on the frequency domain
voltage response data $V(k)$, it is possible to detect whether the
electrochemical system behaves as an LTI, NLTI, LTV or NLTV system, as
illustrated in Fig. 12.

Figure 12: Detection of nonlinear and nonstationary behaviour in measured
data. An odd random phase multisine signal (blue crosses) with a detector for
odd nonlinear distortions is applied for $P=10$ periods at two different
excitation amplitudes. In the LTI case, the voltage response only has spectral
content at DC (dots at $f=0\,\mathrm{Hz}$) and at the excited frequencies. For
a higher excitation amplitude, the system might become NLTI, and spectral
content appears at integer multiples of the fundamental frequency. When the
system is nonstationary, hyperbolic shapes become visible around the excited
frequencies, and in the nonlinear case around all the integer multiples of the
fundamental frequency.

LTI:

$V(k)$ has only spectral content at DC and the excited harmonics
$P\mathbb{H}_{\mathrm{exc}}$.

NLTI:

$V(k)$ has spectral content at DC and the excited frequencies, and also at
harmonics that are integer multiples of the fundamental frequency $1/T_{p}$,
i.e., $P\mathbb{H}_{\mathrm{nl}}$ with
$\mathbb{H}_{\mathrm{nl}}=\\{0,1,2,...\\}$. Both even and odd nonlinearities
are detected, since there is spectral content at the even left out harmonics
($0$, $2$, $4$, $6$) and at the odd left out harmonic ($5$).

LTV:

$V(k)$ consists of hyperbolic-like shapes around DC and the excited
frequencies. These shapes, called _skirts_ in the literature [122], are due to
the smooth time-varying function modulating the multisine components.

NLTV:

$V(k)$ consists of hyperbolic-like shapes around DC, the excited frequencies
and the nonlinear harmonics.

In Fig. 12 noiseless data are considered. However, real-life measurements also
contain noise that is distributed over all DFT frequencies. Still, it is easy
to distinguish between nonlinearities, nonstationarities and noise. This can,
for instance, be seen from Fig. 11 where the measured voltage spectrum clearly
satisfies the LTI constraints and the noise is at least $1000$ times smaller
than the linear response. For a _nonzero-mean_ multisine excitation, the
response will likely be nonstationary. This would look like a sampled version
of Fig. 8.

From all the spectra illustrated in Fig. 12, only the LTI one should be
processed with classical EIS techniques, as discussed in the previous section.
If only nonlinearities are detected, one can estimate the BLA (see Section
8.2). If time-variation is detected, but no nonlinear distortions, time-
varying impedance data can be estimated from the data record (see Section 9),
if there are nonlinear distortions too, but they are limited, the BLTVA can be
estimated from the data using operando EIS (see Section 9.3).

###  7.2 Advantages of an ORP multisine excitation

By measuring the response of the electrochemical system to an ORP multisine
excitation for a large number of periods, we can easily detect the presence of
nonstationarities and nonlinearities, and estimate the noise level. It is
possible to do this over a wide frequency range, in one single experiment, as
studied in [38, 39, 40, 42]. The advantages of using a multisine excitation
over a single-sine one have already been discussed in depth (Section 3.2.3).
However, an interested reader might wonder why it is particularly advantageous
to use a multisine with _random phases_ and only _odd_ excited harmonics. The
reasons are as follows:

The random phases are chosen to minimise the crest factor of the multisine.
Note that there exist optimal ways to achieve this, however, using random
phases has the advantage of simple implementation [57].

The choice to excite only _odd_ harmonics, on the other hand, is more
involved. In fact, the advantage is twofold. First, it allows distinguishing
between even and odd nonlinear distortions. The nonlinear distortions are
present at DFT lines which are the product of the excited harmonic numbers and
the degrees of the relevant monomials in the Volterra series [104]. When
exciting only odd harmonics, for odd nonlinear behaviour (odd degree monomials
in the Volterra series), the nonlinear distortions are present at odd
harmonics ($\mathrm{odd}\times\mathrm{odd}=\mathrm{odd}$), while for even
nonlinear behaviour (even degree monomials in the Volterra series) they are
present at even harmonics ($\mathrm{odd}\times\mathrm{even}=\mathrm{even}$).
If we would only excite even harmonics, it would not be possible to
distinguish between the two ($\mathrm{even}\times\mathrm{odd}=\mathrm{even}$
and $\mathrm{even}\times\mathrm{even}=\mathrm{even}$). When exciting all
harmonics, even and odd, it is also not possible to distinguish between the
types of nonlinearities.

The second advantage is that even nonlinearities do not have a contribution at
the excited frequencies. Hence, when computing the BLA or BLTVA using a
multisine excitation, even nonlinear distortions introduce no uncertainty, as
further discussed in Section 8.2.

##  8 Nonlinear impedance estimation

###  8.1 Leading order nonlinear impedance estimation

It is appropriate to perform NLEIS (Section 4.1.2) in the frequency domain.
For this purpose, we apply a zero mean single-sine excitation, that is,
$i_{0}(t)=0$, $M=1$ and $\mathbb{H}_{\text{exc}}=\\{1\\}$ in (47), measured
over an integer number of periods $P\geq 1$. When choosing the right amplitude
$I_{1}$, as discussed in Section 4.1.2, we obtain estimates of the leading
order nonlinear impedance coefficients as,

$\displaystyle\hat{Z}_{h,h}(\omega_{P})=\frac{V(hP)}{I(P)^{h}}$
$\displaystyle\omega_{P}=\frac{2\pi P}{T}=\frac{2\pi}{T_{p}}.$ (61)

Measuring a higher number of periods $P>1$ introduces an averaging effect of
the stochastic noise in the experiments. Different frequencies $\omega_{P}$
are applied _sequentially_. NLEIS with multisine excitation should also be
possible, however not yet investigated. One has to be careful for leaving
enough gaps in the excitation signal to see the integer harmonics in the
response.

###  8.2 Best linear approximation

Independently of the excitation amplitude of _single-sine_ experiments, the
BLA can be estimated from (38), where different frequencies can be applied
_sequentially_ ,

$\displaystyle\hat{Z}(\omega_{P})=\frac{V(P)}{I(P)}$
$\displaystyle\omega_{P}=\frac{2\pi P}{T}=\frac{2\pi}{T_{p}}.$ (62)

For _multisine_ excitations measured for an integer number of periods $P$, the
BLA impedance estimate is computed exactly as for classical impedance
experiments (59). If an ORP multisine is used, that is, exciting odd
harmonics, only odd nonlinearities will introduce spectral content at the
excited frequencies. Indeed, even nonlinear distortions can only be present at
even multiples of the odd harmonics, resulting in even harmonics. Accordingly,
the BLA will have uncertainties due to noise and also due to odd nonlinear
distortions that introduce spectral content at other odd harmonics (see [26,
Section 3.4, p. 78]). The effect of the noise can again be reduced by
measuring a larger number of periods. The odd nonlinear distortions, on the
other hand, do not depend on the number of measured periods. They only depend
on the system and the excitation, that is, the amplitudes and excited
frequencies. The uncertainty of the BLA can be reduced by exciting fewer
frequencies in the multisine. Recall that the BLA depends on the excitation
(38), however, using Riemann-equivalent multisines, one can adapt the number
of excited frequencies without changing the BLA [123].

##  9 Time-varying impedance estimation

Estimating the time-varying impedance $Z(\omega_{k},t)$ cannot simply be done
by dividing DFT spectra. Different approaches have been developed over the
last two decades to reveal the time-varying impedance from data. Both single-
sine and multisine techniques exist. However, in this review paper we
intentionally restrict to multisine excitations, since we believe this is the
only correct solution (see Section 4.2.3). With ever-growing computation
power, more and more complex techniques have been developed. There are two
main approaches. One implies the use of windowing/filtering of the current and
voltage data in time or frequency domain, respectively, and gives an average
value of the impedance inside the selected time frame. Another uses the
mathematical regression of the voltage response spectrum to extract the time-
variation. These approaches are detailed next, chronologically.

###  9.1 FFT-EIS applied to nonstationary data

The first attempt to obtain time-varying impedance data was by Bond, Schwall
and Smith in 1977 [43, 124, 125]. This simple and intuitive approach is an
extension of Smith’s FFT-EIS discussed in Section 6.

Instead of applying only a zero-mean excitation signal, the excitation signal
is now superimposed on a slower cyclic voltammetry excitation, introducing
nonstationarity in the measured system. While collecting current and voltage
data, the FFT-EIS is performed on short subrecords. This is a type of
windowing. Accordingly, for each subrecord, which corresponds to a certain
point in the voltage trajectory of the cyclic voltammetry, the time-averaged
impedance of the subrecord is computed using the techniques of Section 6.
Basically, classical impedance measurements are applied to nonstationary data,
but in very small time-windows such that the impedance could be assumed
constant within each subrecord.

Since the period length is inversely proportional to the frequency, the
impedance is only measured at high frequencies, such that short subrecords can
be taken. Results of the impedance (or admittance) varying over the cyclic
voltammetry trajectory are only shown for two particular frequencies above
$300\,\mathrm{Hz}$. This is a strong limitation of the technique.

Later, in 1997, FFT-EIS was implemented on a micro computer [126], for a
multisine superimposed onto a staircase DC ramped voltage excitation. Time-
varying impedance measurements could be obtained within the more broadband
frequency range of $50\,\mathrm{Hz}$ – $50\,\mathrm{kHz}$. Later, Sacci and
Harrington [46, 47], developed measurement apparatus to obtain time-varying
impedance data using FFT-EIS with multisine excitation superimposed on a
cyclic voltammogram.

##### Comment

It is noteworthy that obtaining time-varying impedance data is much easier at
high frequencies. This on the one hand due to low frequency noise (so-called
$1/f$ noise), but also due to the drift signal (or also called trend)
$v_{0}(t)$ in (41) having a decreasing shape over frequency, and hence, hiding
low-frequency content [63, 59]. Moreover, for logarithmically distributed
excited frequencies, nonstationarities are more easily detected at the high
frequencies since the excited frequencies are more separated (expressed in DFT
lines), thus making the ‘skirts’ better visible. Remarkably, time-varying
behaviour is mostly present at the low frequencies, this is possibly due to
mass transport and/or charge transfer kinetics. Also, when measuring at high
frequencies, the measurement time is usually shorter, therefore
nonstationarities are smaller.

###  9.2 Time-frequency analysis methods

During the nineties, time-frequency analysis, as described by L. Cohen555Not
to confuse with the great Canadian singer-songwriter. [127], became a widely
used tool in signal processing. Time-frequency analysis describes how the
spectral content of a signal $x(t)$ is changing in time, which is exactly
needed for time-varying impedance estimation. The workhorse for this job is
the short-time Fourier transform (STFT), which computes the Fourier transform
of a signal restricted by a window function $w(t)$,

$\displaystyle\mathrm{STFT}\\{x\\}(\omega,t)$
$\displaystyle=\int_{-\infty}^{\infty}w(t^{\prime}-t)x(t^{\prime})e^{-j\omega
t^{\prime}}\mathrm{d}t^{\prime}$
$\displaystyle=\mathcal{F}\\{w(t^{\prime}-t)x(t^{\prime})\\},$ (63)

with the Fourier transform acting on the variable $t^{\prime}$. The most
commonly used window functions are the Gaussian, Hamming and Blackmann-Harris
windows. These windows reach their largest values in the center, and decrease
smoothly towards the borders.

####  9.2.1 STFT-EIS

After the millenium change, Darowicki [44, 128, 45] proposed to estimate the
time-varying impedance under multisine excitation as the ratio of the STFT of
voltage and current,

$\displaystyle
Z(\omega,t)=\frac{\mathcal{F}\\{w(t^{\prime}-t)v(t^{\prime})\\}}{\mathcal{F}\\{w(t^{\prime}-t)i(t^{\prime})\\}},$
(64)

with again the Fourier transforms acting on the variable $t^{\prime}$.
Assuming the window $w(t)$ to be a symmetric function centered around zero,
the impedance at time $t$ is computed by selecting the time domain data around
this time-instant, and dividing the corresponding spectra of voltage and
current. Note that FFT-EIS on nonstationary data is a special case of STFT-
EIS, with a rectangular window $w(t)$.

In practice, of course, only discrete-time data (50) is available. The
impedance can then be computed by [45],

$\displaystyle\hat{Z}(\omega_{k},t_{n})=\frac{V(k,n)I^{*}(k,n)}{I(k,n)I^{*}(k,n)}$
$\displaystyle t_{n}=nT_{s},$ (65a) with the DFT acting on subrecords of
$N_{w}$ data points centered around $n$, which are windowed by the function
$w(t)$, $\displaystyle X(k,n)$
$\displaystyle=\frac{1}{N_{w}}\sum_{n^{\prime}=n-N_{w}/2}^{n+N_{w}/2-1}w(n^{\prime}-n)x(n^{\prime})e^{-j2\pi
kn^{\prime}/N_{w}},$ (65b)

$x=i,v$. A division of the cross- and auto spectra is chosen here since the
signals are not periodic anymore, and this estimator is recommended for random
excitations (see Section 6). The time-varying impedance estimate is computed
at the harmonics $k$ where the numerator of (65a) shows peak values,
corresponding to the excited frequencies of the multisine. To make these peaks
visible and not overlapping too much, it is important that enough periods are
measured. Moreover, the choice of the window $w(t)$ is crucial in this
technique.

####  9.2.2 Dynamic Multi-Frequency Analysis

Later, Battistel and La Mantia [48, 107, 54] proposed the dynamic multi-
frequency analysis (DMFA) for estimating time-varying impedance data. Here,
the time-varying impedance is computed by filtering the current and voltage
spectra around the excited frequencies of a multisine, and taking the inverse
Fourier transform,

$\displaystyle
Z(\omega,t)=\frac{\mathcal{F}^{-1}\\{W(\omega^{\prime}-\omega)V(\omega^{\prime})\\}}{\mathcal{F}^{-1}\\{W(\omega^{\prime}-\omega)I(\omega^{\prime})\\}},$
(66)

with the inverse Fourier transforms acting on $\omega^{\prime}$. The function
$W(\omega)$ implements a filtering operation. This process is called
_quadrature filtering_ since only the spectrum at positive frequencies is
considered and the inverse Fourier transforms give complex-valued results.

For sampled frequency domain data (54), with multisine excitation, the time-
varying impedance data estimates translate into

$\displaystyle\hat{Z}(\omega_{k},t_{n})=\frac{v(k,n)}{i(k,n)},$ (67a) with for
$k\in P\mathbb{H}_{\text{exc}}$ the inverse DFT acting on frequency domain
subrecords of $N_{W}$ data points centered around $k$, which are filtered by
the function $W(f)$, $\displaystyle
x(k,n)=\sum_{k^{\prime}=k-N_{W}/2}^{k+N_{W}/2-1}W(k^{\prime}-k)X(k^{\prime})e^{j2\pi
k^{\prime}n/N_{W}}\quad x=i,v.$ (67b)

Here also, it is important that enough periods are measured, such that around
the excited frequencies, enough frequency domain data is available to extract
the time-variation.

####  9.2.3 Equivalence between STFT-EIS and DMFA

For a symmetrical window, $w(t)=w(-t)$, and $W(\omega)=\mathcal{F}\\{w(t)\\}$,
the continuous definitions of the STFT-EIS (64) and the DMFA (66) are
equivalent (as proven in E),

$\displaystyle\frac{\mathcal{F}\\{w(t^{\prime}-t)v(t^{\prime})\\}}{\mathcal{F}\\{w(t^{\prime}-t)i(t^{\prime})\\}}=\frac{\mathcal{F}^{-1}\\{W(\omega^{\prime}-\omega)V(\omega^{\prime})\\}}{\mathcal{F}^{-1}\\{W(\omega^{\prime}-\omega)I(\omega^{\prime})\\}}.$
(68)

Since the Fourier and inverse Fourier transforms act on, respectively,
$t^{\prime}$ and $\omega^{\prime}$, both left and right hand side of (68) are
a function of $\omega$ and $t$. Note that both latter equations for extracting
the time-varying impedance are heuristic, and do not exactly match with the
theoretical definition (4.2.1). Nonetheless, these approximations may be
accurate enough in practice.

Since the STFT-EIS and DMFA approaches have a mathematically equivalent
definition of the impedance, the difference boils down to the choice of the
window $w(t)$, or equivalently the filter $W(\omega)$, and the actual
implementation. The properties of the symmetrical window or filter can mainly
be studied by its width. This width can, for instance, be defined by the
variances

$\displaystyle\sigma^{2}_{t}=\frac{\int_{-\infty}^{\infty}t^{2}|w(t)|^{2}\mathrm{d}t}{\int_{-\infty}^{\infty}|w(t)|^{2}\mathrm{d}t}\
\text{and}\
\sigma^{2}_{\omega}=\frac{\int_{-\infty}^{\infty}\omega^{2}|W(\omega)|^{2}\mathrm{d}\omega}{\int_{-\infty}^{\infty}|W(\omega)|^{2}\mathrm{d}\omega}.$
(69)

Similarly as in quantum mechanics where the uncertainty principle prohibits to
measure simultaneously the position and velocity of an electron with arbitrary
precision, also in this case it is not possible to measure the impedance with
arbitrary precision in both time and frequency. Accordingly, the so-called
_Gabor limit_ [127] states that

$\displaystyle\sigma^{2}_{t}\sigma^{2}_{\omega}\geq\frac{1}{4}.$ (70)

The time-selectivity $\sigma^{2}_{t}$ and frequency resolution
$\sigma^{2}_{\omega}$ cannot both be made arbitrarily small. One has to trade-
off between one and the other.

##### STFT-EIS

Darowicki [44] uses a Gaussian window, which has the property that its Fourier
transform is Gaussian as well,

$\displaystyle w(t)=e^{-\frac{\lambda}{2}t^{2}}\ \iff\
W(\omega)=\sqrt{\frac{2\pi}{\lambda}}e^{-\frac{\omega^{2}}{2\lambda}}.$ (71)

For this choice of window, the Gabor limit reduces to

$\displaystyle\sigma^{2}_{t}\sigma^{2}_{\omega}=\frac{1}{4}\quad\text{with}\quad\lambda=\frac{1}{2\sigma^{2}_{t}}=2\sigma^{2}_{\omega}.$
(72)

Accordingly, an increase of the time selectivity leads to a deterioration of
the frequency resolution, and vice versa. The time and frequency resolution
are determined by the hyper parameter $\lambda$. The larger $\lambda$, the
more resolution in time, but the less resolution in frequency, resulting in a
trade-off.

An issue with STFT-EIS is that first subrecords are taken, followed by the
windowing. As a consequence, we do not actually reach the Gabor limit. Also
the DFT is taken on subrecords, which obviously contain less data points than
the total record, hence, as discussed in Section 5, the SNR is poorer. Another
issue is obtaining time-varying impedance data at low frequencies. For being
able to measure low frequencies, the width of the window should, at least,
comprise a period of this frequency. The period length is inversely
proportional to the frequency, hence, for measuring low frequencies, the
window width should be large, deteriorating the time resolution. On the other
hand, since STFT-EIS is applied to time windows, the advantage is that the
time-varying impedance estimation can be done in real-time.

##### DMFA

Battistel and La Mantia [107], on the other hand, directly define the
quadrature filter,

$\displaystyle
W(\omega)=\frac{\big{(}1+e^{-q^{2}}\big{)}^{2}}{\big{(}1+e^{-q\frac{\omega+\Delta\omega}{\Delta\omega}}\big{)}\big{(}1+e^{q\frac{\omega-\Delta\omega}{\Delta\omega}}\big{)}},$
(73)

where $q$ is a factor determining the roll-off of the filter and
$\Delta\omega$ is its bandwidth. Note that in the limits, we have that

$\displaystyle\lim_{q\rightarrow 0}W(\omega)=1$ (74a) and
$\displaystyle\lim_{q\rightarrow\infty}W(\omega)=\left\\{\begin{array}[]{ll}1&\mbox{if
}-\Delta\omega\leq\omega\leq\Delta\omega\\\ 0&\mbox{else.}\end{array}\right.$
(74d)

The objective of this filter is to mimic a rectangular filter, while being
continuous.

The advantage of the DMFA over the STFT-EIS is that only one DFT should be
performed on the entire time-series data record (50) to obtain the frequency
domain data (54). Here, the measurement noise is distributed over all the DFT
lines, resulting in a higher SNR (see Section 5). The time-varying impedance
data is then directly obtained by applying the inverse DFT to small windowed
subrecords around the excited frequencies of the multisine. This has
advantages in computation time (see [107] Section 2.3). Also the width of the
filter can be chosen in function of the spacing of the excited frequencies,
possibly different for each excited frequency. Moreover, an analysis of the
influence of measurement noise on the time-varying impedance data was
performed [129]. However, the time-varying impedance estimation cannot be done
in real-time.

The main advantage of these time-frequency analysis methods is the accessible
implementation. However, both time-frequency analysis methods have difficulty
in estimating time-varying impedance data at low frequencies, this is due to
the drift signal. Moreover, they do not account for treating nonlinear
distortions, even though in the DMFA this could be implemented. These problems
are each solved by operando EIS, as detailed next.

###  9.3 Operando EIS

Operando EIS, as developed by Hallemans et al. [49, 63, 59], is an extension
of ORP-EIS. Here, we estimate the time-varying impedance by using the
definition of the BLTVA (46) as a model structure. Note that if no nonlinear
distortions are present in the measurements, the BLTVA and the time-varying
impedance in (4.2.1) are equal. Nonlinearities in the measurements are
detected, quantified and classified, and the noise level is estimated [49].
Also uncertainty bounds are included on the estimated impedance data.
Moreover, drift signals are suppressed, allowing to access low frequencies
[63].

The idea is to write the time-varying impedance as a truncated series
expansion in a set of known basis functions in time,

$\displaystyle Z(\omega,t)=\sum_{p=0}^{N_{p}}Z_{p}(\omega)b_{p}(t).$ (75)

Figure 13: The first four Legendre polynomials in time (a) and frequency
domain (b).

The basis functions $b_{p}(t)$ are chosen as Legendre polynomials (F and Fig.
13 (a)), since these benefit from good numerical conditioning [122]. Using
(75), the frequency and time dependencies are separated. Since the basis
functions are known, only the impedances $Z_{p}(\omega)$ should be estimated.

The time-varying nonlinear distortions $v_{\mathrm{s}}(t)$ are also expanded
in series,

$\displaystyle
v_{\mathrm{s}}(t)=\sum_{p=0}^{N_{p}}v_{\mathrm{s},p}(t)b_{p}(t),$ (76)

with $v_{\mathrm{s},p}(t)$ the time-invariant nonlinear distortions generated
by a Volterra series (23), meaning that
$V_{\mathrm{s},p}(\omega)=\mathcal{F}\\{v_{\mathrm{s},p}(t)\\}$ is only
nonzero at the integer multiples of the fundamental frequency of the periodic
excitation. Plugging (75) and (76) in the excitation-response relation of NLTV
systems (46) yields,

$\displaystyle
v(t)=v_{0}(t)+\sum_{p=0}^{N_{p}}\mathcal{F}^{-1}\\{Z_{p}(\omega)I_{\mathrm{exc}}(\omega)+V_{\mathrm{s},p}(\omega)\\}b_{p}(t).$
(77)

The drift signal $v_{0}(t)$ is unknown and hides low frequency content.
Therefore, it should also be modelled, for instance by Legendre polynomials
[42, 49],

$\displaystyle v_{0}(t)=\sum_{q=0}^{N_{q}}\theta_{q}b_{q}(t).$ (78)

Drift signals can also be removed by differencing, as detailed in [63]. This
has better performance, however, the mathematics are more involved, so for
this review paper we restrict to modelling the drift signal with basis
functions.

Taking the DFT of (77) gives,

$\displaystyle
V(k)=V_{0}(k)+\sum_{p=0}^{N_{p}}\Big{(}Z_{p}(\omega_{k})I_{\mathrm{exc}}(k)+V_{\mathrm{s},p}(k)\Big{)}\ast
B_{p}(k),$ (79)

with

$\displaystyle V_{0}(k)=\sum_{q=0}^{N_{q}}\theta_{q}B_{q}(k).$ (80)

Here it was used that a product in the time domain becomes a convolution in
the frequency domain and $B_{p}(k)$ is the DFT of the Legendre polynomials,
shown in Fig. 13 (b). For a multisine excitation measured over an integer
number of periods $P$, we have that $I_{\mathrm{exc}}(k)$ is only nonzero at
the harmonics $P\mathbb{H}_{\mathrm{exc}}$ and $V_{\mathrm{s},p}(k)$ is only
nonzero at the harmonics $P\mathbb{H}_{\mathrm{nl}}$, with
$\mathbb{H}_{\mathrm{nl}}=\\{0,1,2,3,...\\}$. Note that
$\mathbb{H}_{\mathrm{exc}}\subset\mathbb{H}_{\mathrm{nl}}$ . Accordingly, (79)
can be simplified as,

$\displaystyle
V(k)=\sum_{q=0}^{N_{q}}\theta_{q}B_{q}(k)+\sum_{p=0}^{N_{p}}\sum_{k^{\prime}\in
P\mathbb{H}_{\mathrm{nl}}}\theta_{p}(k^{\prime})B_{p}(k-k^{\prime}),$ (81)

with

$\displaystyle\theta_{p}(k^{\prime})=Z_{p}(\omega_{k^{\prime}})I_{\mathrm{exc}}(k^{\prime})+V_{\mathrm{s},p}(k^{\prime}).$
(82)

Eq. (81) is linear in the parameters $\theta$, hence, the data can be written
in matrix form,

$\displaystyle V=K\theta+N_{v},$ (83)

where $V$ is a stacked vector of the measured voltage spectra $V(k)$, the
regression matrix $K$ consists of regressors $B_{p}(k)$ centered around the
harmonics $P\mathbb{H}_{\mathrm{nl}}$, the parameter vector $\theta$ contains
the parameters $\theta_{q}$, $q=0,1,...,N_{q}$, $\theta_{p}(k^{\prime})$,
$p=0,1,...,N_{p}$ and $k^{\prime}\in P\mathbb{H}_{\mathrm{nl}}$, and $N_{V}$ a
vector representing the noise. The optimal parameters are estimated in linear
least squares sense,

$\displaystyle\hat{\theta}$
$\displaystyle=\arg\min_{\theta}(V-K\theta)^{H}(V-K\theta)$ (84a)
$\displaystyle=(K^{H}K)^{-1}K^{H}V.$ (84b)

For long data records, the regression problem becomes too large to solve in
one go. Therefore, it is proposed to solve it in local frequency bands as
detailed in [49]. We retrieve the impedance and nonlinear distortion estimates
from the estimated parameter vector $\hat{\theta}$,

$\displaystyle\hat{Z}_{p}(\omega_{k})$
$\displaystyle=\frac{\hat{\theta}_{p}(k)}{I_{\mathrm{exc}}(k)}$ $\displaystyle
k\in P\mathbb{H}_{\mathrm{exc}}$ (85a)
$\displaystyle\hat{V}_{\mathrm{s},p}(k)$ $\displaystyle=\hat{\theta}_{p}(k)$
$\displaystyle k\in
P\big{(}\mathbb{H}_{\mathrm{nl}}\setminus\mathbb{H}_{\mathrm{exc}}\big{)}$
(85b)

Finally, the estimate of the time-varying impedance at the excited frequencies
is obtained as,

$\displaystyle\hat{Z}(\omega_{k},t)=\sum_{p=0}^{N_{p}}\hat{Z}_{p}(\omega_{k})b_{p}(t)$
$\displaystyle k\in P\mathbb{H}_{\mathrm{exc}}.$ (86)

Odd nonlinear distortions and noise introduce uncertainties on the time-
varying impedance estimates $\hat{Z}(\omega_{k},t)$. For the computation of
uncertainty bounds due to nonlinear distortions and noise, the reader is
referred to [49], where also a noise estimation is performed.

The strength of operando EIS is that it computes the BLTVA of NLTV data,
together with uncertainty bounds. Hence, electrochemical systems not
satisfying linearity can still be monitored using an impedance. However, the
odd nonlinearities must not be too strong. Moreover, since the drift signal is
modelled as well, low frequency information becomes available, which is
important for some applications. Therefore, it is applicable to a wide range
of experiments. Also, it does measure the actual definition of the time-
varying impedance, or the BLTVA in the nonlinear case.

The trade-off between time- and frequency resolution, however, remains. To be
able to extract the time-variation, we should leave enough empty DFT lines in
between the excited lines (corresponding to measuring a large integer number
of periods), which decreases the resolution of the excited lines.

##  10 A case study on commercial Li-ion batteries

Li-ion batteries are chosen as a case study to illustrate the important
concepts in this article on real-life measurements. Similar experiments are
performed as in [59]. Measurements are performed on a pristine Samsung
INR21700-48X cell, placed in a thermal chamber at $5^{\circ}$C or
$25^{\circ}$C. The commercially available Samsung 48X is a
$4.8\text{\,}\mathrm{A}\mathrm{h}$ $21$ $700$ cell format with cathodes based
on lithiated metal oxide (Co, Ni, Al) and anodes based on intercalation
graphite and blended Si.

Current and voltage data are collected using the Gamry Interface 5000E
potentiostat. Besides running classical EIS experiments, this potentiostat
allows to apply user-defined excitations, and measure current and voltage
data. The sampling frequency for these user-defined excitations is limited to
$200\text{\,}\mathrm{H}\mathrm{z}$. The current range that can be applied is
limited to $[-5,5]$ A.

An odd random phase multisine signal $i_{\mathrm{exc}}(t)$ is designed with
period $T_{p}=3$ min. The $76$ excited frequencies are chosen as odd
harmonics, $\mathbb{H}_{\text{exc}}=\\{1,3,5,...\\}$, logarithmically
distributed between $f_{\mathrm{min}}=5.6\,\mathrm{mHz}$ and
$f_{\mathrm{max}}=80\,\mathrm{Hz}$. The phases are chosen randomly, uniformly
distributed in $[0,2\pi)$. Since noise is often more prominent at low
frequencies, and while doing operando experiments the low frequency content is
hidden by drift signals, the amplitudes of the multisine are chosen with a
decreasing shape over frequency, with root-mean-square (RMS) of $0.8$ A rms.

###  10.1 Estimating classical impedance data

For _classical_ EIS experiments, the battery is first entirely discharged,
then charged at a $C/3$ rate (constant current of $4.8/3=1.6$ A) until the
desired SOC level ($10,20,\ldots,90$%). Two hours of relaxation are allowed
such that the battery reaches steady-state and the voltage reaches the OCV
value. Then, the zero-mean excitation signal $i(t)=i_{\mathrm{exc}}(t)$ is
looped for $P=10$ periods, that is for $T=PT_{p}=30$ min. The measured voltage
at the different operating points ($\text{SOC}=10,30,50,70,90$% and
$\mathrm{T}=25^{\circ}$C) is shown in Fig. 14. Note that indeed, the voltage
data looks periodic with period $3$ min, and that the data is nicely centered
around the OCV value, as should be the case in LTI measurements. Since the
measurements are periodic, the classical impedances are easily computed from
Section 6 (59).

The spectra of the measured current and voltage, and estimated impedance for
the $10$% SOC and $25^{\circ}$C operating point are shown in Fig. 11. Note
that for the chosen excitation amplitudes, the battery behaves very linearly.
No nonlinear distortions nor nonstationarities can be detected. However, a
(small) drift signal is present. Also, noise is clearly present in the
measurements, but it is fairly low, an SNR of at least $1000$ is obtained.

The estimated impedances at different operating points (depending on SOC and
temperature) are shown in Fig. 3. We do indeed obtain a different impedance
for each of the operating points. Note that if we want to perform quicker
experiments, we can measure less periods, for instance $P=4$, leading to
experiments of $T=12$ min, however, with a higher noise floor.

Figure 14: Voltage response of the Samsung 48X cell in LTI conditions at
$25^{\circ}$C. The zero-mean excitation signal $i_{\mathrm{exc}}(t)$ is used
and $P=10$ periods of the current and voltage are measured at different SOC
levels. Every color indicates a measurement at a different SOC. Note that for
every SOC, the OCV is also different. The time-series are subsampled $150$
times.

###  10.2 Estimating time-varying impedance data

For obtaining _time-varying_ impedance data, the battery is charged with a
$C/2$ current, with the multisine superimposed,

$\displaystyle i(t)=2.4\,\mathrm{A}+i_{\mathrm{exc}}(t).$ (87)

The top graph of Fig. 15 shows the measured current and voltage data at
$5^{\circ}$C. Due to the DC offset of $2.4$ A in the current, the battery is
charging, and the voltage goes up, leading to a drift signal superimposed on
the multisine response. The measurement is stopped when the voltage hits the
safety bound of $4.2$ V. For a constant current charging of $C/2$ we would
expect the $4.2$ V to be reached after $2$ h. However, due to the multisine
added on top of this charging current the safety limit is reached prematurely.
Accordingly, $P=29$ and $P=31$ periods of the excitation could be measured for
the $5^{\circ}$C and $25^{\circ}$C experiments, respectively.

The middle graph of Fig. 15 shows the SOC, with values from $0$ % to $72.5$ %,
and the battery’s surface temperature, which increases slightly due to the
charging current. This was also shown in the SOC-temperature plane in Fig. 9.

The spectra of the current and voltage data of the $5^{\circ}$C experiment are
shown in Fig. 16. The bottom graph shows the entire spectra with a logarithmic
frequency axis, while the top graph shows zoomed spectra in different
frequency bands, each $36\,\mathrm{mHz}$ wide, with a linear frequency axis.
Note the general decreasing shape of the drift spectrum $V_{0}(k)$ in the
voltage spectra which hides the low frequency content. For the lowest zoomed
frequency band (top left), the time-invariant contributions at the excited
frequencies barely exceed the drift spectrum, and the skirts are completely
hidden. At frequencies close to $1\,\mathrm{Hz}$, the skirt around the excited
frequency is a little more visible, still the drift spectrum hides
information. At frequencies close to $80\,\mathrm{Hz}$ the skirts are clearly
visible.

The time-varying impedance at $1\,\mathrm{Hz}$, estimated using operando EIS
[59], is shown in the bottom graph of Fig. 15. The time-varying impedance at
all excited frequencies is shown in Fig. 9. Even though the drift signal hides
the low frequency content, clean impedance data can be obtained at these
frequencies using operando EIS. Note that the impedance is highest at low SOC,
and that the impedance while charging is different from while resting, as also
observed in [59].

Figure 15: Experiment performed on a Samsung 48X cell in time-varying
conditions in a thermal chamber at $5^{\circ}$C. Top graph: current excitation
and voltage response. The current has a DC offset of $2.4$ A, hence, the
battery charges, and the voltage increases. Middle graph: SOC, obtained by
Coulomb counting, and the external parameter temperature during the
experiment. Since the battery is charging with a constant current plus zero-
mean multisine, the SOC increases linearly and the temperature increases
slightly. Bottom graph: time-variation of the impedance at
$0.9389\,\mathrm{Hz}$ obtained from operando EIS [59]. Figure 16: Current and
voltage spectra of an experiment performed on a Samsung 48X cell in time-
varying conditions at $5^{\circ}$C. Top graph: spectra of current and voltage
in three different frequency bands of each $38.5\,\mathrm{mHz}$ wide, with a
linear frequency axis. Note the decreasing shape of the drift spectrum hiding
low frequency content. Bottom graph: entire spectra, with a logarithmic
frequency axis.

###  10.3 Estimating equivalent circuit model parameters

Studying the processes going on at different time-scales in the impedance data
is often performed by mapping the data onto ECM parameters. For the measured
battery impedance data on Samsung 48X cells, we choose the ECM of Fig. 17.
This ECM can be linked to the SPM for batteries [34]. The resistance $R_{0}$
(yellow) is related to the electrolyte resistance, the first $RC$-branch with
Warburg element (purple) is related to the diffusion and the second
$RC$-branch (green) to the electrochemical kinetics. The corresponding
parametric impedance yields,

$\displaystyle
Z_{\mathrm{ECM}}(\omega,\theta)=R_{0}+Z_{C_{1}}\text{//}(R_{1}+Z_{\mathrm{W}})+Z_{C_{\mathrm{ct}}}\text{//}R_{\mathrm{ct}},$
(88a) where the parameter vector $\theta$ is given by
$\displaystyle\theta=[R_{0},R_{1},C_{1},R_{\mathrm{ct}},C_{\mathrm{ct}},\mathrm{W},\alpha],$
(88b) the symbol ‘//’ stands for the parallel connection, that is,
$\displaystyle
Z_{X}(\omega)\text{//}Z_{Y}(\omega)=\frac{Z_{X}(\omega)Z_{Y}(\omega)}{Z_{X}(\omega)+Z_{Y}(\omega)},$
(88c) and the impedance of a capacitor and Warburg element, respectively,
yield, $\displaystyle Z_{C}(\omega)=\frac{1}{Cj\omega}\quad\text{and}\quad
Z_{\mathrm{W}}(\omega)=\frac{\mathrm{W}}{(j\omega)^{\alpha}}.$ (88d)

The Nyquist chart in Fig. 17 illustrates the contribution of the three
branches in series (yellow, purple and green) on the total impedance (black),
being the sum of the three other colors.

The ECM parameters $\theta$ can now be estimated from impedance data by
minimising the cost function,

$\displaystyle\hat{\theta}=\arg\min_{\theta}\sum_{k\in
P\mathbb{H}_{\text{exc}}}|\hat{Z}(\omega_{k})-Z_{\mathrm{ECM}}(\omega_{k},\theta)|^{2}.$
(89)

This cost function is nonlinear in the parameters $\theta$, hence, a nonlinear
solver is required. Here, we use a hybrid of _particle swarm optimisation_
[130, 131, 132] and the built-in MATLAB function lsqnonlin. Fits over the
frequency band [16.7 mHz,50 Hz] are obtained with mean relative errors over
frequency all smaller than $0.2$% and $0.5$% for the classical and time-
varying impedance data, respectively.

The estimated ECM parameters for the measured impedance data of the Samsung
48X cell at $5^{\circ}$C and $25^{\circ}$C are shown in Fig. 17. The ECM
parameters of the classical impedance data at different operating points are
shown as dots, while the ECM parameters of the time-varying impedance data
along trajectories are shown as continuous lines. The temperature is assumed
approximately constant during the experiments. It is observed that the
parameters obtained in operating and classical conditions are not necessarily
equal. For Li-ion batteries, this was already studied in [56], where the
charge transfer resistance $R_{\mathrm{ct}}$ at a certain SOC is smaller while
charging than while resting.

As an application, Zhu et al. [58] propose a fast charging protocol by
applying a charging current inversely proportional to the time-varying charge-
transfer resistance $R_{\mathrm{ct}}$, tracked using operando EIS.

Equivalent circuit model

Corresponding Nyquist plot and estimated parameters

Figure 17: Equivalent circuit model and estimated parameters for the Samsung
48X cell. The Nyquist plot shows the influence of the different branches
(yellow, purple and green) on the total impedance (black). In the ECM
parameter graphs, dots represent the parameters obtained from classical EIS at
different operating points, while the continuous lines show the parameters
obtained from time-varying impedance data along operating trajectories. For
the time-varying experiments, the temperature is assumed approximately
constant at $5$ or $25^{\circ}$C.

##  11 Conclusions & outlook

Classical EIS provides impedance data of electrochemical systems at selected
frequencies. Due to the constraints of linearity and stationarity, the
impedance data is only valid for small amplitude excitations and at fixed
operating points. Nonetheless, measuring classical impedance data is a
powerful tool for monitoring electrochemical systems.

Models beyond linearity and stationarity, such as nonlinear leading-order and
time-varying impedance, reveal higher-dimensional impedance data, valid over
larger excitation amplitudes and along operating trajectories. This higher-
dimensional impedance data contains additional information to classical
impedance data, which is promising for electrochemical applications. One
could, for instance, increase the accuracy of health forecasting of Li-ion
batteries using nonlinear and/or time-varying impedance data as indicator.
This is also extendable to other electrochemical applications, such as
detecting corrosion or studying coatings.

It is shown that the multisine excitation is a strong asset for modelling
electrochemical systems. It allows nonlinear and nonstationary behaviour to be
detected from the measured current and voltage data. If this current and
voltage data does not satisfy the linearity and stationarity constraints,
higher-dimensional impedance data can still be extracted.

## Acknowledgements

This article is dedicated to Rik Pintelon. NH has received funding from the
Eutopia mobility programme. NH and JL are supported financially by the Fund
for Scientific Research (FWO Vlaanderen) and the Flemish government, Belgium
(grant number: METH1). FLM has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement n° 772579). This work was supported by the
Fraunhofer Internal Programs under Grant No. Attract 028-602604 ProLIBs.

##  Appendix A The Fourier transform

$\displaystyle X(\omega)$ $\displaystyle=\mathcal{F}\\{x(t)\\}$
$\displaystyle=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}\mathrm{d}t$ (90)
$\displaystyle x(t)$ $\displaystyle=\mathcal{F}^{-1}\\{X(\omega)\\}$
$\displaystyle=\int_{-\infty}^{\infty}X(\omega)e^{j\omega t}\mathrm{d}\omega.$
(91)

##  Appendix B Volterra series coefficients

Second order nonlinearity

$\displaystyle V_{2,0}$
$\displaystyle=\frac{1}{4}\big{(}Z_{2}(-\omega,\omega)+Z_{2}(\omega,-\omega)\big{)}I^{2}$
(92) $\displaystyle V_{2,1}$ $\displaystyle=0$ (93) $\displaystyle V_{2,2}$
$\displaystyle=\frac{1}{2}Z_{2}(\omega,\omega)I^{2}$ (94)

Third order nonlinearity

$\displaystyle V_{3,0}$ $\displaystyle=0$ (95) $\displaystyle V_{3,1}$
$\displaystyle=\frac{1}{8}\big{(}Z_{3}(-\omega,\omega,\omega)+Z_{3}(\omega,-\omega,\omega)+Z_{3}(\omega,\omega,-\omega)\big{)}I^{3}$
(96) $\displaystyle V_{3,2}$ $\displaystyle=0$ (97) $\displaystyle V_{3,3}$
$\displaystyle=\frac{1}{8}Z_{3}(\omega,\omega,\omega)I^{3}$ (98)

##  Appendix C Linearising an NLTI system around an operating trajectory

The origin of the nonstationarity could be proven by applying a particular
excitation to a Volterra series, consisting of a slow part, dictating the
trajectory, and a fast part, which is the excitation,

$\displaystyle
i(t)=\underbrace{i_{0}(t)}_{\text{slow}}+\underbrace{i_{\text{exc}}(t)}_{\text{fast}}.$
(99)

As an example, the slow part could be a positive constant current for charging
a battery, and the fast part a multisine. By assuming that the fast
perturbation has a small amplitude, and hence, only the linear part of the
Volterra series ($n=1$) is needed with respect to $i_{\text{exc}}(t)$, the
voltage response $v(t)$ can also be separated into a slow and fast part,

$\displaystyle v(t)=v_{0}(t)+v_{\text{exc}}(t),$ (100a) with $\displaystyle
v_{0}(t)$
$\displaystyle=\mathrm{OCV}+\sum_{n=1}^{n_{\mathrm{max}}}\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}z_{n}(\tau_{1},...,\tau_{n})\prod_{l=1}^{n}i_{0}(t-\tau_{l})\mathrm{d}\tau_{l}$
$\displaystyle v_{\text{exc}}(t)$
$\displaystyle=\int_{-\infty}^{\infty}\underbrace{z(\tau,t)}_{\text{depends on
}z_{n}\text{'s and }i_{0}(t)}i_{\text{exc}}(\tau)\mathrm{d}\tau.$ (100b)

The slow part $v_{0}(t)$ is called the drift signal, and solely depends on the
slow excitation $i_{0}(t)$. The fast response is now the convolution of a two-
dimensional impulse response $z(\tau,t)$ with the excitation. This two-
dimensional impulse response function explicitly depends on the time of
excitation $t$, such that stationarity is not satisfied anymore. Moreover,
this function is shown to depend on the generalised impulse responses
$z_{n}(\tau_{1},\dots,\tau_{n})$ and the slow signal $i_{0}(t)$.

##  Appendix D The discrete Fourier transform

$\displaystyle X(k)$
$\displaystyle=\frac{1}{N}\sum_{n=0}^{N-1}x(nT_{s})e^{-j\frac{2\pi kn}{N}}$
(101) $\displaystyle x(nT_{s})$
$\displaystyle=\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi kn}{N}}$ (102)

##  Appendix E Equivalence between (64) and (66)

Using the following properties of the Fourier transform,

$\displaystyle X(\omega)=\mathcal{F}\\{x(t)\\}\qquad x=x,y,w,i,v$ (103a)
$\displaystyle\mathcal{F}\\{x(t)y(t)\\}=X(\omega)\ast
Y(\omega)=\int_{-\infty}^{\infty}X(\omega-\omega^{\prime})Y(\omega^{\prime})\mathrm{d}\omega^{\prime}$
(103b) $\displaystyle\mathcal{F}\\{x(t^{\prime}-t)\\}=X(\omega)e^{-j\omega
t},$ (103c)

where $\mathcal{F}$ acts on $t^{\prime}$, one finds that,

$\displaystyle\mathcal{F}\\{w(t^{\prime}-t)x(t^{\prime})\\}$
$\displaystyle=\int_{-\infty}^{\infty}W(\omega-\omega^{\prime})e^{-j(\omega-\omega^{\prime})t}X(\omega^{\prime})\mathrm{d}\omega^{\prime}$
$\displaystyle=e^{-j\omega
t}\int_{-\infty}^{\infty}W(\omega-\omega^{\prime})X(\omega^{\prime})e^{j\omega^{\prime}t}\mathrm{d}\omega^{\prime}$
$\displaystyle=e^{-j\omega
t}\mathcal{F}^{-1}\\{W(\omega-\omega^{\prime})X(\omega^{\prime})\\}.$ (104)

Assuming that $w(t)=w(-t)$, one has that $W(\omega)=W(-\omega)$, accordingly,

$\displaystyle\frac{\mathcal{F}\\{w(t^{\prime}-t)v(t^{\prime})\\}}{\mathcal{F}\\{w(t^{\prime}-t)i(t^{\prime})\\}}$
$\displaystyle=\frac{\mathcal{F}^{-1}\\{W(\omega-\omega^{\prime})V(\omega^{\prime})\\}}{\mathcal{F}^{-1}\\{W(\omega-\omega^{\prime})I(\omega^{\prime})\\}}$
$\displaystyle=\frac{\mathcal{F}^{-1}\\{W(\omega^{\prime}-\omega)V(\omega^{\prime})\\}}{\mathcal{F}^{-1}\\{W(\omega^{\prime}-\omega)I(\omega^{\prime})\\}}$
(105)

with all Fourier and inverse Fourier transforms acting on, respectively,
$t^{\prime}$ and $\omega^{\prime}$.

##  Appendix F Legendre polynomials

The Legendre polynomials $L_{p}(x)$, $p=0,1,...$, $x\in[-1,1]$ are the
solution of Legendre’s differential equation

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}x}\Big{(}(1-x^{2})\frac{\mathrm{d}L_{p}(x)}{\mathrm{d}x}\Big{)}+p(p+1)L_{p}(x)=0.$
(106)

The basis functions $b_{p}(t)$ are chosen as rescaled Legendre polynomials
over the interval $[0,T]$, that is,

$\displaystyle b_{p}(t)=L_{p}\Big{(}\frac{2t}{T}-1\Big{)}.$ (107)

## References

  * [1] M. E. Orazem, B. Tribollet, Electrochemical Impedance Spectroscopy, Wiley, 2008 (2008).
  * [2] J. R. MacDonald, Impedance Spectroscopy, Wiley, 1987 (1987).
  * [3] C. Gabrielli, Identification of electrochemical processes by frequency response analysis, Technical Report Number 04/83 2 (1984) 25 (1984).
  * [4] C. Gabrielli, Once upon a time there was EIS, Electrochimica Acta 331 (2020) 135324 (2020).
  * [5] S. Wang, J. Zhang, O. Gharbi, V. Vivier, M. Gao, M. E. Orazem, Electrochemical impedance spectroscopy, Nature Reviews Methods Primers 1 (1) (2021) 1–21 (2021).
  * [6] M. Doyle, T. F. Fuller, J. Newman, Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell, Journal of the Electrochemical society 140 (6) (1993) 1526 (1993).
  * [7] C.-H. Chen, F. B. Planella, K. O’regan, D. Gastol, W. D. Widanage, E. Kendrick, Development of experimental techniques for parameterization of multi-scale lithium-ion battery models, Journal of The Electrochemical Society 167 (8) (2020) 080534 (2020).
  * [8] F. B. Planella, M. Sheikh, W. D. Widanage, Systematic derivation and validation of a reduced thermal-electrochemical model for lithium-ion batteries using asymptotic methods, Electrochimica Acta 388 (2021) 138524 (2021).
  * [9] F. Mansfeld, Electrochemical impedance spectroscopy (EIS) as a new tool for investigating methods of corrosion protection, Electrochimica Acta 35 (10) (1990) 1533–1544 (1990).
  * [10] F. Mansfeld, Use of electrochemical impedance spectroscopy for the study of corrosion protection by polymer coatings, Journal of applied electrochemistry 25 (3) (1995) 187–202 (1995).
  * [11] R. I. Revilla, B. Wouters, F. Andreatta, A. Lanzutti, L. Fedrizzi, I. De Graeve, EIS comparative study and critical equivalent electrical circuit (EEC) analysis of the native oxide layer of additive manufactured and wrought 316l stainless steel, Corrosion Science 167 (2020) 108480 (2020).
  * [12] U. Tröltzsch, O. Kanoun, H.-R. Tränkler, Characterizing aging effects of lithium ion batteries by impedance spectroscopy, Electrochimica acta 51 (8-9) (2006) 1664–1672 (2006).
  * [13] D. Andre, M. Meiler, K. Steiner, C. Wimmer, T. Soczka-Guth, D. Sauer, Characterization of high-power lithium-ion batteries by electrochemical impedance spectroscopy. i. experimental investigation, Journal of Power Sources 196 (12) (2011) 5334–5341 (2011).
  * [14] W. Waag, S. Käbitz, D. U. Sauer, Experimental investigation of the lithium-ion battery impedance characteristic at various conditions and aging states and its influence on the application, Applied energy 102 (2013) 885–897 (2013).
  * [15] R. R. Richardson, P. T. Ireland, D. A. Howey, Battery internal temperature estimation by combined impedance and surface temperature measurement, Journal of Power Sources 265 (2014) 254–261 (2014).
  * [16] C. Pastor-Fernández, K. Uddin, G. H. Chouchelamane, W. D. Widanage, J. Marco, A comparison between electrochemical impedance spectroscopy and incremental capacity-differential voltage as Li-ion diagnostic techniques to identify and quantify the effects of degradation modes within battery management systems, Journal of Power Sources 360 (2017) 301–318 (2017).
  * [17] N. Meddings, M. Heinrich, F. Overney, J.-S. Lee, V. Ruiz, E. Napolitano, S. Seitz, G. Hinds, R. Raccichini, M. Gaberšček, et al., Application of electrochemical impedance spectroscopy to commercial Li-ion cells: A review, Journal of Power Sources 480 (2020) 228742 (2020).
  * [18] X. Zhu, L. F. Macía, J. Jaguemont, J. de Hoog, A. Nikolian, N. Omar, A. Hubin, Electrochemical impedance study of commercial $\mathrm{LiNi_{0.80}Co_{0.15}Al_{0.05}O_{2}}$ electrodes as a function of state of charge and aging, Electrochimica Acta 287 (2018) 10–20 (2018).
  * [19] Y. Zhang, Q. Tang, Y. Zhang, J. Wang, U. Stimming, A. A. Lee, Identifying degradation patterns of lithium ion batteries from impedance spectroscopy using machine learning, Nature communications 11 (1) (2020) 1–6 (2020).
  * [20] M. Gaberšček, Understanding Li-based battery materials via electrochemical impedance spectroscopy, Nature Communications 12 (1) (2021) 1–4 (2021).
  * [21] P. K. Jones, U. Stimming, A. A. Lee, Impedance-based forecasting of lithium-ion battery performance amid uneven usage, Nature Communications 13 (1) (2022) 1–9 (2022).
  * [22] Z. He, F. Mansfeld, Exploring the use of electrochemical impedance spectroscopy (EIS) in microbial fuel cell studies, Energy & Environmental Science 2 (2) (2009) 215–219 (2009).
  * [23] S. M. R. Niya, M. Hoorfar, Study of proton exchange membrane fuel cells using electrochemical impedance spectroscopy technique–a review, Journal of Power Sources 240 (2013) 281–293 (2013).
  * [24] T. Kailath, Linear systems, Prentice-Hall Englewood Cliffs, NJ, 1980 (1980).
  * [25] L. Ljung, System Identification: Theory for the User, Prentice Hall, 1999 (1999).
  * [26] R. Pintelon, J. Schoukens, System Identification: a Frequency Domain Approach, John Wiley & Sons, 2012 (2012).
  * [27] M. Schetzen, The Volterra and Wiener theories of nonlinear systems (1980).
  * [28] S. Boyd, L. Chua, Fading memory and the problem of approximating nonlinear operators with Volterra series, IEEE Transactions on circuits and systems 32 (11) (1985) 1150–1161 (1985).
  * [29] N. Harting, N. Wolff, F. Röder, U. Krewer, Nonlinear frequency response analysis (NFRA) of lithium-ion batteries, Electrochimica Acta 248 (2017) 133–139 (2017).
  * [30] F. Fasmin, R. Srinivasan, Nonlinear electrochemical impedance spectroscopy, Journal of The Electrochemical Society 164 (7) (2017) H443 (2017).
  * [31] M. D. Murbach, V. W. Hu, D. T. Schwartz, Nonlinear electrochemical impedance spectroscopy of lithium-ion batteries: experimental approach, analysis, and initial findings, Journal of The Electrochemical Society 165 (11) (2018) A2758 (2018).
  * [32] T. Vidaković-Koch, T. Miličić, L. A. Živković, H. S. Chan, U. Krewer, M. Petkovska, Nonlinear frequency response analysis: a recent review and perspectives, Current Opinion in Electrochemistry 30 (2021) 100851 (2021).
  * [33] M. A. Zabara, G. Katırcı, B. Ülgüt, Non-linear harmonics in EIS of batteries with lithium anodes: Proper controls and analysis, Electrochimica Acta (2022) 140969 (2022).
  * [34] T. L. Kirk, A. Lewis-Douglas, D. Howey, C. P. Please, S. J. Chapman, Nonlinear electrochemical impedance spectroscopy for lithium-ion battery model parameterization, Journal of The Electrochemical Society (2022).