#  Nonlinear electrochemical impedance spectroscopy for lithium-ion battery
model parameterization

Toby L. Kirk  Adam Lewis-Douglas  David Howey  Colin P. Please  S. Jon Chapman

###### Abstract

In this work we analyse the local nonlinear electrochemical impedance
spectroscopy (NLEIS) response of a lithium-ion battery and estimate model
parameters from measured NLEIS data. The analysis assumes a single-particle
model including nonlinear diffusion of lithium within the electrode particles
and asymmetric charge transfer kinetics at their surface. Based on this model
and assuming a moderately-small excitation amplitude, we systematically derive
analytical formulae for the impedances up to the second harmonic response,
allowing the meaningful interpretation of each contribution in terms of
physical processes and nonlinearities in the model. The implications of this
for parameterization are explored, including structural identifiability
analysis and parameter estimation using maximum likelihood, with both
synthetic and experimentally measured impedance data. Accurate fits to
impedance data are possible, however inconsistencies in the fitted diffusion
timescales suggest that a nonlinear diffusion model may not be appropriate for
the cells considered. Model validation is also demonstrated by predicting
time-domain voltage response using the parameterized model and this is shown
to have excellent agreement with measured voltage time-series data
($11.1\text{\,}\mathrm{m}\mathrm{V}$ RMSE).

###### keywords:

battery , model , impedance , parameter , nonlinear

††journal: Journal of the Electrochemical Society

[inst1]organization=Mathematical Institute, University of
Oxford,addressline=Andrew Wiles Building, Woodstock Road, city=Oxford,
postcode=OX2 6GG, country=UK

[inst2]organization=The Faraday Institution,addressline=Quad One, Becquerel
Avenue, Harwell Campus, city=Didcot, postcode=OX11 0RA, country=UK
[inst3]organization=Battery Intelligence Lab, Department of Engineering
Science, University of Oxford,city=Oxford, postcode=OX1 3PJ, country=UK

##  1 Introduction

Lithium-ion batteries have emerged as the dominant energy storage solution for
portable electronics and electric vehicles, enabled by their high energy and
power density and decreasing costs [1]. They are also increasingly of interest
for renewable energy integration within the power grid. In all of these
applications, battery models are used to improve understanding of performance
trade-offs at design stage and for operational monitoring of key factors such
as state of charge, temperature, and lifetime. Equivalent circuit models are
widely used in battery management systems, but electrochemical models, such as
those derived from porous electrode theory within the Doyle-Fuller-Newman
(DFN) framework [2], are also increasingly of interest for high-fidelity
predictions in academia and industry.

The effectiveness of electrochemical battery models as a tool for performance
understanding and improvement depends on the realism of their internal
structure and parameters for the particular questions being asked. Even if
measured input-output quantities such as voltage can be predicted accurately,
the internal states may differ significantly from reality. This is compounded
by the large number of parameters required in these models (e.g. the DFN model
has more than 30 parameters) and hence values are often taken from literature,
but may not apply to the specific cell at hand. Fast and noninvasive
characterisation techniques that can completely parameterize a physical model
for a specific cell are therefore desirable.

Electrochemical impedance spectroscopy (EIS) is a widely used noninvasive
frequency-domain characterization technique where a small amplitude sinusoidal
current (or voltage) is applied to a device and the corresponding voltage (or
current) response is measured [3]. The system response is assumed to be
linear, and is traditionally fitted using ‘local’ equivalent circuit models
containing resistors, capacitors, and sometimes other _ad hoc_ components
(e.g. constant phase elements), although this may lead to overfitting and
difficulties in interpretation of the underlying processes—see Ciucci et al.
[4] for a review of the application of EIS to lithium-ion batteries. The
assumption of a linear response in standard EIS necessarily discards
information about nonlinearities that may be useful for model selection and
parameterization. This information may be measured by applying a larger
amplitude sinusoidal excitation to generate harmonics in the response of the
device under test.

The injection of a sufficiently large amplitude sinusoidal current into a
battery will result in a nonlinear voltage response, with harmonic components
not only at the fundamental excitation frequency $\omega$, but at all integer
multiples of this i.e. $0,\omega,2\omega,3\omega,\ldots$; see Fig. 1. The
components $n\omega$ with $n\geq 2$ are referred to as higher harmonics, and
are due entirely to nonlinearities. Several ways of analysing these harmonics
have been considered [5, 6], but one approach is to assume the input amplitude
is small, but larger than for linear EIS. Then a “weakly” nonlinear response
is induced, with analytical predictions from the model possible via a
perturbation expansion. We will refer to this technique as nonlinear EIS
(NLEIS) [7]. It was first widely applied in the context of solid-oxide fuel
cells [8, 9, 10, 11, 12, 13], proving useful for model selection, i.e.
identifying inconsistencies between proposed models and experiment. It has
since been considered for many other electrochemical systems, as discussed in
the review by Fasmin & Srinivasan [14].

Figure 1: Diagram of the NLEIS process, where a sinusoidal current input and
voltage output are transformed to the frequency domain using Fourier
transforms. In practice experimental measurements are pre-processed to remove
transients and windowed (e.g. with a Hann window) before transforming.

There has been limited application of NLEIS to lithium-ion batteries. Murbach
& Schwartz [5], using a DFN-type physical model, calculated the second and
third harmonics by an expansion in current amplitude, but this required
numerical solution of the spatial system of ODEs due to the model complexity.
They showed that the higher harmonics are sensitive to the symmetry of the
charge transfer reaction and to several parameters of interest, e.g. solid-
state and electrolyte diffusivities. The authors also showed qualitative
similarity between measured and predicted nonlinear impedances and that the
second harmonic may be sensitive to cell ageing [15].

There have been several other studies where moderate or large amplitude
sinusoidal currents have been applied to lithium-ion batteries [6, 16, 17,
18]. In particular, nonlinear frequency response analysis (NFRA) is a term
used to refer to a particular type of NLEIS where there is no restriction on
the input amplitude, and only the modulus (magnitude) of the harmonics is
retained. No restriction on input amplitude means a stronger nonlinear
response, however, the model analysis is purely numerical, with direct
simulation in the time domain necessary before transformation to the frequency
domain, which is computationally expensive. The aforementioned “weakly
nonlinear” analysis, on the other hand, can be done directly in the frequency
domain. In addition, the practice in NFRA of discarding the harmonics’ phase
and hence potentially half of the nonlinear information does not seem
justified. NFRA on lithium-ion batteries was explored experimentally by
Harting et al. [6, 16], and signatures of degradation processes (e.g. lithium
plating) were able to be identified. Wolff et al. [17] considered NFRA applied
to a DFN-type model and performed a parameter sensitivity analysis, and they
[18] investigated the signatures of several models of SEI. Large amplitude
multisine signals have recently been considered by Fan et al. [19, 20] who
also performed a sensitivity analysis on a DFN model, giving support to the
notion of asymmetry of charge transfer kinetics reported by others [15].

Previously, NLEIS analysis of battery models has been applied only to
relatively complex physical models with many parameters, such as the DFN
model, making the evaluation and interpretation of the model impedances
difficult. In this paper, we consider instead the NLEIS response of a reduced-
order physical model, the single-particle model (SPM), which is the simplest
physical model that considers each electrode separately and captures the key
battery dynamics and overpotentials, e.g. Butler–Volmer reaction kinetics,
solid-state diffusion of lithium, and double-layer capacitance. This allows us
to derive analytical formulae for the impedances up to the second harmonic,
and interpret each contribution in terms of physical processes, extending the
established language of EIS Nyquist plots to include nonlinear harmonics. The
nonlinearities are identified explicitly, originating from: (i) asymmetric
Butler–Volmer kinetics; (ii) concentration dependence of the exchange current
density; (iii) open circuit potential functions; (iv) solid-state diffusion.
We then consider the implications of this new nonlinear information for
parameter estimation, extending the structural identifiability analysis of a
single-particle model by Bizeray et al. [21], which was limited to linear
impedances. This is complemented with parameter estimates from synthetic and
from experimentally measured NLEIS data using a commercial pouch cell. Our
results show that the model fits the impedance data accurately, capturing many
of the key features of each harmonic and determining all model parameters
excluding open circuit potentials. Lastly, the model is independently
validated in the time domain, showing improvements in accuracy of voltage
predictions over conventional methods.

##  2 Model formulation

###  2.1 The Single-Particle Model (SPM)

The model of a lithium-ion cell we consider is the single-particle model (SPM)
[22, 23, 21]. Here we summarise the model assumptions and governing equations.
The transport of lithium ions in the electrolyte is assumed to be fast, so
that they remain at a constant uniform concentration $c_{e}^{*}$. Throughout,
asterisks denote dimensional quantities. In each electrode, the active
particles that make up the porous electrode medium are all assumed identical
in size and shape. Thus they behave identically, and only a model of a single
representative particle is necessary. A schematic of the model geometry is
shown in Fig. 2.

Figure 2: Single particle model schematic with a representative particle for
each electrode. Not depicted is a resistor in series with the two electrodes.

Let the subscript $\pm$ denote that a quantity applies to, or is defined in,
the positive and negative electrode, respectively. Mass transport of lithium
within spherical particles of radius $\mathcal{R}_{\pm}^{*}$ is modelled here
by spherically-symmetric nonlinear diffusion,

$\displaystyle\frac{\partial c_{\pm}^{*}}{\partial t^{*}}$
$\displaystyle=-\frac{1}{r_{\pm}^{*2}}\frac{\partial}{\partial
r_{\pm}^{*}}\left(r_{\pm}^{*2}N_{r}^{*}\right),\qquad\text{for
}0<r_{\pm}^{*}<\mathcal{R}_{\pm}^{*},$ (1) $\displaystyle N_{r}^{*}$
$\displaystyle=-D_{\pm}^{*}(c_{\pm}^{*})\frac{\partial c_{\pm}^{*}}{\partial
r_{\pm}^{*}},$ (2)

where $c_{\pm}^{*}(r_{\pm}^{*},t^{*})$ is the concentration of lithium in the
electrode material, and $D_{\pm}^{*}(c_{\pm}^{*})$ is the nonlinear solid-
state diffusion coefficient that is a function of concentration. The boundary
conditions are regularity at the particle centres,

$\frac{\partial c_{\pm}^{*}}{\partial r_{\pm}^{*}}=0,\qquad\text{at
}r_{\pm}^{*}=0,$ (3)

and a surface flux, i.e., lithium (de)intercalation into (out of) the
particles,

$-D_{\pm}^{*}(c_{\pm}^{*})\frac{\partial c_{\pm}^{*}}{\partial
r_{\mathrm{\pm}}^{*}}=j_{\pm}^{*},\qquad\text{at
}r_{\mathrm{\pm}}^{*}=\mathcal{R}_{\mathrm{\pm}}^{*},$ (4)

which is modelled by nonlinear Butler–Volmer kinetics,

$\displaystyle
j_{\pm}^{*}=\frac{i_{0,\pm}^{*}(c_{\mathrm{\pm}}^{*})}{2F^{*}}\left[\exp\left(\frac{F^{*}(1-\beta_{\pm})}{R_{g}^{*}T^{*}}\eta_{\mathrm{\pm}}^{*}\right)-\exp\left(-\frac{F^{*}\beta_{\pm}}{R_{g}^{*}T^{*}}\eta_{\mathrm{\pm}}^{*}\right)\right]\qquad\text{at
}r_{\mathrm{\pm}}^{*}=\mathcal{R}_{\mathrm{\pm}}^{*},$ (5) $\displaystyle
i_{0,\pm}^{*}(c_{\mathrm{\pm}}^{*})=m_{\mathrm{\pm}}^{*}(c_{\mathrm{\pm}}^{*})^{\beta_{\pm}}(c_{e}^{*})^{1-\beta_{\pm}}(c_{\pm,\text{max}}^{*}-c_{\pm}^{*})^{1-\beta_{\pm}}.$
(6)

Here $F^{*}$ is Faraday’s constant, $R_{g}^{*}$ is the universal gas constant,
$T^{*}$ is the temperature (assumed constant), $i_{0,\pm}^{*}$ is the exchange
current density (with a dependence on the surface value of
$c_{\mathrm{\pm}}^{*}$, with evaluation on $r_{\pm}^{*}=\mathcal{R}_{\pm}^{*}$
assumed), $m_{\pm}^{*}$ is a reaction rate coefficient,
$c_{\pm,\text{max}}^{*}$ is the maximum lithium concentration in the
positive/negative electrode material. The cathodic transfer coefficient
$\beta_{\pm}$ (and corresponding anodic transfer coefficient $1-\beta_{\pm}$)
is not assumed to be the typical value of $1/2$, but may take any value
between 0 and 1, i.e., $\beta_{\pm}\in(0,1)$. Electrical conductivity in the
electrode materials is assumed large, meaning the electric potential in each
is uniform, only a function of time $t^{*}$, and expressed relative to the
potential in the electrolyte,
$V_{\pm}^{*}(t^{*})=\phi_{\pm}^{*}(t^{*})-\phi_{e}^{*}$. Then
$\eta_{\pm}^{*}(t^{*})=V_{\pm}^{*}(t^{*})-U_{\pm}^{*}(c_{\pm}^{*}|_{r_{\pm}^{*}=R_{\pm}^{*}})$
is the reaction overpotential, where the open circuit potential (OCP) is given
by $U_{\pm}^{*}$, a function of the surface concentration. The evaluation of
$c_{\pm}^{*}$ at $r_{\pm}^{*}=\mathcal{R}_{\pm}^{*}$ in the current density
and overpotential expressions will, for brevity, be understood from here on
and not stated explicitly.

The nonconstant diffusion coefficient $D_{\pm}^{*}(c_{\pm}^{*})$ will be left
as a general function of $c_{\pm}^{*}$ in our analysis, but for our results we
will interpret it as originating from a diffusive flux driven by chemical
potential gradients, such as in concentrated solution theory. Therefore we
take the radial lithium flux to be [24]

$\displaystyle
N_{r}^{*}=\frac{{D}_{\mathrm{Li},\pm}^{*}F^{*}c_{\pm}^{*}}{R_{g}^{*}T^{*}}\frac{\partial
U_{\pm}^{*}}{\partial r_{\pm}^{*}},$ (7)

which has been shown to give more accurate behaviour than a Fickian diffusion
model when parameterized using galvanostatic intermittent titration technique
(GITT). Here ${D}_{\mathrm{Li},\pm}^{*}$ is the constant diffusivity of the Li
species in each electrode material, and the flux (7) gives rise to (1) with
the concentration-dependent diffusivity

$\displaystyle D_{\pm}^{*}(c_{\pm}^{*})$
$\displaystyle=-\frac{{D}_{\mathrm{Li},\pm}^{*}F^{*}}{R_{g}^{*}T^{*}}\frac{\mathrm{d}U_{\pm}^{*}}{\mathrm{d}c_{\pm}^{*}}c_{\pm}^{*}.$
(8)

Note that, as per this equation, the dependence on $c_{\pm}^{*}$ is completely
determined, and depends on the gradient of $U_{\pm}^{*}$.

Finally, conservation of charge relates both of the surface lithium flux
densities to the macroscopic applied circuit current $I_{\text{app}}^{*}$. Let
$a_{\pm}^{*}$ be the surface area per unit volume, related to the electrode
volume fraction $\epsilon_{\pm}$ for spherical particles by
$a_{\pm}^{*}=3\epsilon_{\pm}/\mathcal{R}_{\pm}^{*}$, then the interfacial
current density is related to the applied current via the differential
equation

$C_{\mathrm{dl},\pm}^{*}\frac{\mathrm{d}V_{\pm}^{*}}{\mathrm{d}t^{*}}=\pm\frac{I_{\text{app}}^{*}}{\mathcal{A}^{*}L_{\pm}^{*}a_{\pm}^{*}}-F^{*}j_{\pm}^{*},$
(9)

where $I_{\text{app}}^{*}>0$ corresponds to a charging current,
$C_{\mathrm{dl},\pm}^{*}$ is the double-layer capacitance at the surface of
each electrode, $L_{\pm}^{*}$ is the through-cell thickness of each electrode,
and $\mathcal{A}^{*}$ is the area of the electrodes in the transverse
directions.

We include double-layer capacitance here due to its relevance in modelling the
frequency domain behaviour, in particular at high frequencies. However, a
model neglecting double-layer capacitance has computational advantages and is
particularly useful in time-domain simulations of scenarios where short
timescales are less important, e.g., charging or discharging. We may neglect
capacitance by setting $C_{\mathrm{dl},\pm}^{*}=0$ in (9), reducing it to the
algebraic equation

$j_{\pm}^{*}=\pm\frac{I_{\text{app}}^{*}}{F^{*}\mathcal{A}^{*}L_{\pm}^{*}a_{\pm}^{*}}.$
(10)

If the applied current is prescribed, then substituting (10) into the boundary
condition (4) allows (1)-(4) to be solved for $c_{\mathrm{\pm}}^{*}$, and (5)
subsequently inverted to give the potentials in each electrode (relative to
the electrolyte potential)

$V_{\pm}^{*}=U_{\pm}^{*}(c_{\pm}^{*})+\eta_{\pm}^{*},$ (11)

With the potentials $V_{\pm}^{*}$ given either by the solution of (9) (when
capacitance is included) or the explicit expression (11) (when capacitance is
neglected), the terminal voltage for the cell is then

$\displaystyle V^{*}(t^{*})$
$\displaystyle=\phi_{\mathrm{+}}^{*}-\phi_{-}^{*}+R_{s}^{*}I_{\text{app}}^{*}=V_{+}^{*}-V_{-}^{*}+R_{s}^{*}I_{\text{app}}^{*}.$
(12)

A series resistance, $R_{s}^{*}$, is included to account for linear
resistances arising from, e.g., current collector contacts, and ionic and
electronic conduction in the electrodes and electrolyte.

###  2.2 Scaling and nondimensional SPM

To reduce the number of parameters required and establish the minimal set of
physical parameter groups necessary to use the model, we proceed to
nondimensionalise it. The internal variables $c_{\pm}^{*}$ and $r_{\pm}^{*}$
are scaled by the maximum lithium concentrations $c_{\pm,\mathrm{max}}^{*}$
and particle radii $\mathcal{R}_{\pm}^{*}$, respectively. The external
variables, i.e., those directly measurable via the cell terminals such as the
current $I_{\mathrm{app}}^{*}$, voltage $V^{*}$, and time $t^{*}$, are scaled
with typical values $J^{*}$, $\Phi^{*}$, $\tau^{*}$. It is convenient for
these values to be known _a priori_ , that is, not chosen from some internal
property of the cell itself. We choose the potential scale to be the thermal
voltage, $\Phi^{*}=F^{*}/(R_{g}^{*}T^{*})$, and $J^{*}=1$ A, $\tau^{*}=1$ s.

The nondimensional variables are defined as

$\displaystyle c_{\pm}$ $\displaystyle=c_{\pm}^{*}/c_{\pm,\mathrm{max}}^{*},$
$\displaystyle r$ $\displaystyle=r_{\pm}^{*}/\mathcal{R}_{\pm}^{*},$
$\displaystyle t$ $\displaystyle=t^{*}/\tau^{*},$ $\displaystyle
D_{\pm}(c_{\pm})$
$\displaystyle=\frac{\tau^{*}D_{\pm}^{*}(c_{\pm}^{*})}{\mathcal{R}_{\pm}^{*2}}$
(13) $\displaystyle I$ $\displaystyle=I_{\mathrm{app}}^{*}/J^{*},$
$\displaystyle(V,V_{\pm},U_{\pm},\eta_{\pm})$
$\displaystyle=(V^{*},V_{\pm}^{*},U_{\pm}^{*},\eta_{\pm}^{*})/\Phi^{*},$
$\displaystyle j_{\pm}$
$\displaystyle=\frac{\tau^{*}j_{\pm}^{*}}{c_{\pm,\mathrm{max}}^{*}\mathcal{R}_{\pm}^{*}},$
(14)

resulting in the nondimensional model equations

$\displaystyle\frac{\partial c_{\pm}}{\partial t}$
$\displaystyle=\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(D_{\pm}(c_{\pm})r^{2}\frac{\partial c_{\pm}}{\partial r}\right),$
$\displaystyle\text{for }0<r<1,$ (15) $\displaystyle\frac{\partial
c_{\pm}}{\partial r}$ $\displaystyle=0,$ $\displaystyle\text{at }r=0,$ (16)
$\displaystyle-D_{\pm}(c_{\pm})\frac{\partial c_{\pm}}{\partial
r_{\mathrm{\pm}}}$ $\displaystyle=j_{\pm},$ $\displaystyle\text{at }r=1,$ (17)
$\displaystyle c_{\pm}$ $\displaystyle=c_{\pm,\mathrm{init}}$
$\displaystyle\text{at }t=0,$ (18)

and

$C_{\pm}\frac{\mathrm{d}V_{\pm}}{\mathrm{d}t}=\pm
I(t)-\frac{j_{\pm}}{\xi_{\pm}},$ (19)

with surface fluxes, overpotentials, and terminal voltage given by

$\displaystyle j_{\pm}$
$\displaystyle=\frac{2\xi_{\pm}}{R_{\pm}(c_{\pm})}E_{\pm}(\eta_{\pm}),$ (20)
$\displaystyle\eta_{\pm}$ $\displaystyle=V_{\pm}-U_{\pm}(c_{\pm}),$ (21)
$\displaystyle V$ $\displaystyle=V_{+}-V_{-}+R_{s}I(t),$ (22)

and $C_{\pm}$ defined later in (28). For convenience we have introduced the
notation

$\displaystyle E_{\pm}(\eta)$
$\displaystyle=\frac{1}{2}\left[\exp\left((1-\beta_{\pm})\eta\right)-\exp\left(-\beta_{\pm}\eta\right)\right],$
(23) $\displaystyle R_{\pm}(c)$
$\displaystyle=\frac{2\chi_{\pm}}{c^{\beta_{\pm}}(1-c)^{1-\beta_{\pm}}},$ (24)

where $\chi_{\pm}$ is defined below in (28). If the reactions are symmetric,
i.e. $\beta_{\pm}=1/2$, then $E_{\pm}(\eta)$ reduces to $\sinh(\eta/2)$. Here
$R_{\pm}$ is the reciprocal of the (nondimensional) exchange current density,
and can be interpreted as a nonlinear reaction resistance, with a dependence
on the surface concentration $c_{\pm}$ (evaluation at
$r_{\pm}=\mathcal{R}_{\pm}$ understood). This dependence has the nonlinear,
but assumed known, functional form in (24) with parameters $\beta_{\pm}$ and
$\chi_{\pm}$. The factor of 1/2 in (20) is necessary for $R_{\pm}$ to
correspond to the well-known charge transfer resistance if the model is
linearized—see the Nonlinear Impedance section.

If we neglect double-layer capacitance, then we may set $C_{\pm}=0$ in (19),
which reduces to

$j_{\pm}=\pm\xi_{\pm}I(t),$ (25)

allowing $j_{\pm}$ to be eliminated from the model. In this case, substituting
for $j_{\pm}$ in (20) and rearranging for the overpotential gives

$\displaystyle\eta_{\pm}$
$\displaystyle=E_{\pm}^{-1}\left[\pm\frac{1}{2}R_{\pm}(c_{\pm})I(t)\right],$
(26)

where $E_{\pm}^{-1}$ denotes the inverse of $E_{\pm}$, which has no convenient
closed form representation in general (note that if $\beta_{\pm}=1/2$ then
$E_{\pm}^{-1}(x)=2\sinh^{-1}x$). The electrode potentials $V_{\pm}$ and
terminal voltage $V$ are then calculated from (21)-(22).

The remaining nondimensional parameter groups are

$\displaystyle D_{\pm}(c_{\pm})$
$\displaystyle=-\frac{\tau^{*}}{\tau_{\mathrm{d},\pm}^{*}}\frac{\mathrm{d}U_{\pm}}{\mathrm{d}c_{\pm}}c_{\pm},$
$\displaystyle\xi_{\pm}$
$\displaystyle=\frac{J^{*}\tau^{*}}{3Q_{\mathrm{th},\pm}^{*}},$
$\displaystyle\chi_{\pm}$
$\displaystyle=\frac{R_{\mathrm{ct,typ},\pm}^{*}J^{*}}{\Phi^{*}},$ (27)
$\displaystyle C_{\pm}$
$\displaystyle=\frac{C_{\mathrm{dl},\pm}^{*}\Phi^{*}\mathcal{A}^{*}L_{\pm}^{*}a_{\pm}^{*}}{J^{*}\tau^{*}},$
$\displaystyle R_{s}$ $\displaystyle=\frac{R_{s}^{*}J^{*}}{\Phi^{*}},$ (28)

which are defined in terms of the following physically meaningful quantities
of each electrode

$\displaystyle\tau_{\mathrm{d},\pm}^{*}$
$\displaystyle=\frac{\mathcal{R}_{\pm}^{*2}}{{D}_{\mathrm{Li},\pm}^{*}}=\text{solid-
state diffusion timescales}$ (29) $\displaystyle R_{\mathrm{ct,typ,\pm}}^{*}$
$\displaystyle=\frac{\Phi^{*}}{\mathcal{A}^{*}L_{\pm}^{*}a_{\pm}^{*}m_{\pm}^{*}(c_{e}^{*})^{1-\beta_{\pm}}c_{\pm,\mathrm{max}}^{*}}=\text{typical
charge transfer resistances}$ (30) $\displaystyle Q_{\mathrm{th},\pm}^{*}$
$\displaystyle=F^{*}\epsilon_{\pm}c_{\pm,\mathrm{max}}^{*}L_{\pm}^{*}\mathcal{A}^{*}=\text{theoretical
electrode capacities}$ (31)

To recognise that $R_{\mathrm{ct,typ,\pm}}^{*}$ is an order-of-magnitude
estimate of charge transfer resistance, notice that
$\mathcal{A}^{*}L_{\pm}^{*}a_{\pm}^{*}$ is the total active surface area in
each electrode, and
$m_{\pm}^{*}(c_{e}^{*})^{1-\beta_{\pm}}c_{\pm,\mathrm{max}}^{*}$ is the
typical exchange current density, i.e. magnitude of $i_{0,\pm}^{*}$ from 6.
Thus, the denominator in (30) is the typical reaction current.

###  2.3 Model parameters

Here we summarize the model parameters necessary to use the nondimensional SPM
(15)-(24).

####  2.3.1 OCPs and electrode balancing

These are quantities that relate to the rest-state of the cell, and should be
determined by investigations at equilibrium. They are:

  * 1.

$U_{\pm}(c_{\pm})$, the electrode OCPs (relative to a Li/Li+ reference) as a
function of stoichiometry $c_{\pm}\in[0,1]$;

  * 2.

$\xi_{\pm}$, the ratio of the theoretical electrode capacities
$Q_{\mathrm{th},\pm}^{*}$ to the reference scale $J^{*}\tau^{*}$;

  * 3.

$c_{\pm,\mathrm{init}}$, the initial electrode stoichiometries. These can be
expressed in terms of the depth-of-discharge (DoD) of the cell, or
$Q=Q^{*}/Q_{\mathrm{cap}}^{*}\in[0,1]$, where $Q^{*}$ is the discharge
capacity and $Q_{\mathrm{cap}}^{*}$ is the rated cell capacity. Then
$c_{\pm,\mathrm{init}}$ can be calculated using the linear relations

$\displaystyle c_{\pm}$
$\displaystyle=c_{\pm}^{0\%}+(c_{\pm}^{100\%}-c_{\pm}^{0\%})Q,$ (32)

where $c_{\pm}^{0\%}$ and $c_{\pm}^{100\%}$ are the stoichiometries at 0% and
100% DoD (i.e. $Q=0$ and $Q=1$). Note that these can be related to the
theoretical electrode capacities, and hence $\xi_{\pm}$, via

$\displaystyle c_{\pm}^{100\%}-c_{\pm}^{0\%}$
$\displaystyle=\pm\frac{Q_{\mathrm{cap}}^{*}}{Q_{\mathrm{th},\pm}^{*}}=\pm\frac{3Q_{\mathrm{cap}}^{*}}{J^{*}\tau^{*}}\xi_{\pm},$
(33)

therefore only one of the stoichiometries in each electrode, say
$c_{\pm}^{0\%}$, must be provided independently.

In general, this amounts to requiring the knowledge of 2 functions
($U_{\pm}(c_{\pm})$) and 4 constants ($\xi_{\pm}$, $c_{\pm}^{0\%}$) related to
the balancing of the electrodes. These are typically found invasively, by
disassembling a similar cell and reassembling into two half cells and then
measuring the OCVs—, or by the use of reference electrodes. For the
experimental data used in this study, the $U_{\pm}$ are taken to be known as a
function of DoD (or $Q$), measured using a minimally invasive reference
electrode [21, 25]—see the Methods section. We note, however, that the use of
a reference electrode does not allow the electrode OCPs to be determined
outside the normal operation range of the full cell. To do so, half-cell
measurements would be required.

####  2.3.2 Dynamical parameters

With the electrode balancing parameters determined, the remaining 9 parameter
groups necessary for the dynamic modelling of the cell (nonzero currents)
consist of 4 per electrode and 1 series resistance, as follows:

  * 1.

$\tau_{\mathrm{d},\pm}=\tau_{\mathrm{d},\pm}^{*}/\tau^{*}$, typical diffusion
timescale (scaled by the reference timescale $\tau^{*}=1$ s);

  * 2.

$\chi_{\pm}$, typical nondimensional charge transfer resistances;

  * 3.

$\beta_{\pm}$, cathodic charge transfer coefficients;

  * 4.

$C_{\pm}$, nondimensional double-layer capacitances;

  * 5.

$R_{s}$, nondimensional series resistance.

This list of parameters can be written as a vector,

$\displaystyle\bm{\theta}$
$\displaystyle=(\tau_{\mathrm{d},+}\,\,\,\chi_{+}\,\,\,\beta_{+}\,\,\,C_{+}\,\,\,\tau_{\mathrm{d},-}\,\,\,\chi_{-}\,\,\,\beta_{-}\,\,\,C_{-}\,\,\,R_{s})^{T}.$
(34)

Structural identifiability and estimation of these parameters from NLEIS data
will be explored as part of the results of this paper, in the Results section.

###  2.4 Model nonlinearities

Since NLEIS is a nonlinear extension of EIS, it is useful to list the various
sources of nonlinearity present in the model, namely

  * 1.

the OCPs, i.e. the nonlinear functions $U_{\pm}(c_{\pm})$ which encode much of
the thermodynamic behaviour of the electrodes, and

  * 2.

the Butler–Volmer reaction kinetics. These contribute nonlinearities in two
different ways: (i) the exponential dependence of the reaction current on
overpotential, given by $E_{\pm}(\eta_{\pm})$ in (23); (ii) the dependence of
the exchange current density, or equivalently the charge transfer resistance
$R_{\pm}(c_{\pm})$ in (24), on lithium concentration. Finally,

  * 3.

nonlinear diffusion in the particles, i.e. the dependence of the diffusion
coefficient on lithium concentration.

Each of these nonlinearities contributes to the measured nonlinear impedance,
and in the following section we will derive explicit formulae for the
impedance, making each distinct contribution apparent.

##  3 Nonlinear impedance

###  3.1 General description of NLEIS

The technique of (galvanostatic) NLEIS involves the application of a
sinusoidal current, $I(t;\omega,\hat{I})=\hat{I}\mathrm{e}^{\mathrm{i}\omega
t}+\hat{I}\mathrm{e}^{-\mathrm{i}\omega
t}=\text{Re}\\{2\hat{I}\mathrm{e}^{\mathrm{i}\omega t}\\}$ to a device under
test, where, by choice of phase, $\hat{I}$ is real and positive, giving the
real representation $I=2\hat{I}\cos\omega t$ (it is more convenient to work
with $\hat{I}$, rather than the real amplitude $2\hat{I}$). Then, the
nonlinear voltage response of the system, $V(t;\omega,\hat{I})$, is measured.
After the decay of any initial transients, $V$ will in general be periodic in
$t$ with period $P=2\pi/\omega$ matching the input current, and hence have the
Fourier representation

$\displaystyle V(t;\omega,\hat{I})$
$\displaystyle=\sum_{n=-\infty}^{\infty}\widehat{V}_{n}(\omega,\hat{I})\mathrm{e}^{\mathrm{i}n\omega
t},\qquad\text{where }\widehat{V}_{-n}=\overline{\widehat{V}_{n}},$ (35)

and the bar denotes complex conjugation. The component at frequency $\omega$,
matching that of the input current, is referred to as the fundamental.
However, due to nonlinearities in the system, all higher harmonics, i.e. at
frequencies $2\omega,3\omega,\ldots$, are in general also present. The Fourier
coefficients $\widehat{V}_{n}(\omega,\hat{I})$ depend on the excitation
amplitude $\hat{I}$. Thus, we express each Fourier coefficient as a Taylor
series expansion for small current amplitude $\hat{I}\ll 1$, and it can be
shown [15, 14] that in general,

$\widehat{V}_{n}(\omega,\hat{I})=\sum_{r=0}^{\infty}\hat{I}^{n+2r}Z_{n}^{(n+2r)}(\omega)=\hat{I}^{n}Z_{n}^{(n)}+\hat{I}^{n+2}Z_{n}^{(n+2)}+O(\hat{I}^{n+4}),$
(36)

where $Z_{n}^{(n+2r)}(\omega)$, $n,r=0,1,2,\ldots$ are independent of
$\hat{I}$. The subscript refers to the frequency mode or harmonic in the
Fourier series expansion, and the superscript is the index in the Taylor
series expansion in powers of $\hat{I}$. That is, a superscript $(m)$ denotes
a coefficient of $\hat{I}^{m}$. Then $Z_{0}^{(0)}$ is simply the rest (open
circuit) voltage, $Z_{1}^{(1)}$ is the usual linear impedance, and all other
$Z_{n}^{(n+2r)}$ with $n\geq 2$ or $r\geq 1$ are _nonlinear_ impedances.

For a given harmonic $n$ in the voltage, $\hat{I}^{n}Z_{n}^{(n)}$ is the
leading-order response, and so $Z_{n}^{(n)}$ is the leading-order impedance
for that harmonic and can be isolated from $\widehat{V}_{n}$ via

$Z_{n}^{(n)}(\omega)=\lim_{\hat{I}\to
0}\frac{\widehat{V}_{n}(\omega,\hat{I})}{\hat{I}^{n}}.$ (37)

The higher order impedances $Z_{n}^{(n+2r)}$, $r\geq 1$, can be isolated
sequentially in a similar way, provided all lower-order terms for that
harmonic are first subtracted from $\widehat{V}_{n}$. We will focus on the
leading-order impedances in this paper, e.g. $Z_{1}^{(1)}$ and $Z_{2}^{(2)}$,
however $Z_{0}^{(2)}$ is also easily derived in this analysis.

The expansion in current amplitude $\hat{I}\ll 1$ can be performed
analytically on the SPM, solving for each $Z_{n}^{(n+2r)}$, order-by-order, in
a systematic fashion. Formally, all variables
$X\in\\{c_{\pm},\eta_{\pm},j_{\pm},V_{\pm}\\}$ have a series expansion
identical to (35)-(36). Up to $O(\hat{I}^{2})$, which is the highest order we
will consider here, these take the form

$X(t;\omega,\hat{I})=\sum_{n=-\infty}^{\infty}\widehat{X}_{n}(\omega,\hat{I})\mathrm{e}^{\mathrm{i}n\omega
t}\qquad\text{where }\widehat{X}_{-n}=\overline{\widehat{X}_{n}},$ (38)

with Fourier coefficients $\widehat{X}_{n}$, as follows:

$\displaystyle\text{Zero frequency }(0):$ $\displaystyle\widehat{X}_{0}$
$\displaystyle=\widehat{X}_{0}^{(0)}+\hat{I}^{2}\widehat{X}_{0}^{(2)}+O(\hat{I}^{4}),$
(39) $\displaystyle\text{Fundamental/1st harmonic }(\omega):$
$\displaystyle\widehat{X}_{1}$
$\displaystyle=\hat{I}\widehat{X}_{1}^{(1)}+O(\hat{I}^{3}),$ (40)
$\displaystyle\text{2nd harmonic }(2\omega):$ $\displaystyle\widehat{X}_{2}$
$\displaystyle=\hat{I}^{2}\widehat{X}_{2}^{(2)}+O(\hat{I}^{4}),$ (41)
$\displaystyle\text{Higher harmonics }(n\omega,n\geq 3):$
$\displaystyle\vdots$ (42)

###  3.2 SPM neglecting double-layer capacitance

First, we consider the SPM with double-layer capacitance neglected, i.e.,
equations (15)-(18), (23)-(26). In this case there is only a single path that
the current can take through the cell, and the impedance formulae are useful
in expressing the full formulae with capacitance included. We will denote the
impedances without capacitance by a lowercase $z$ to distinguish them from the
full impedances with capacitance included, denoted by $Z$, which will be
stated afterwards.

The details of the expansion are given in section S1 of the supplementary
material, and the resulting impedances for the full cell model can be written
in terms of half-cell impedances $z_{n,\pm}^{(m)}=\widehat{V}_{n,\pm}^{(m)}$,
found from the expansion of the half-cell voltages $V_{\pm}$.

At leading order in the current amplitude, there is only a zero frequency
response, $z_{0,\pm}^{(0)}=U_{\pm}^{(0)}=U_{\pm}(c_{\pm}^{(0)})$, the OCP at
the stoichiometry $c_{\pm}^{(0)}$ for that state of charge.

At first order, $O(\hat{I})$, there is a response at fundamental frequency
given by111The $\pm$ prefactor is due to the fact that a positive (charging)
cell current corresponding to a charging current in the positive electrode but
_discharging_ current in the negative one.

$\displaystyle z_{1,\pm}^{(1)}$
$\displaystyle=\pm\left[R_{\pm}^{(0)}+Z^{\mathrm{W}}_{1,\pm}(\omega)\right],$
(43)

with two contributions: charge transfer resistance $R_{\pm}^{(0)}$ and so-
called finite Warburg impedance due to the diffusion of lithium within a
spherical particle,

$\displaystyle Z^{\mathrm{W}}_{1,\pm}(\omega)$
$\displaystyle=U_{\pm}^{(0)\prime}\frac{\xi_{\pm}}{D_{\pm}^{(0)}}H_{1}\left(\frac{\omega}{D_{\pm}^{(0)}}\right),$
(44)

where a prime denotes a derivative with respect to $c_{\pm}$,

$U_{\pm}^{(0)\prime}=\left.\frac{\mathrm{d}U_{\pm}}{\mathrm{d}c_{\pm}}\right|_{c_{\pm}^{(0)}},$
(45)

and a superscript (0) denotes the leading order (zero current) quantity at
$c_{\pm}^{(0)}$. Note that
$D_{\pm}^{(0)}=-U_{\pm}^{(0)\prime}c_{\pm}^{(0)}/\tau_{\mathrm{d},\pm}$, but
for more generality and brevity we will leave formulae in terms of
$D_{\pm}^{(0)}$. The function $H_{1}(\omega)$ encodes the frequency dependence
of the transport process at this order, and is given by

$\displaystyle H_{1}(\omega)$
$\displaystyle=\frac{\tanh\sqrt{\mathrm{i}\omega}}{\tanh\sqrt{\mathrm{i}\omega}-\sqrt{\mathrm{i}\omega}},$
(46)

with low and high $\omega$ asymptotic behaviours (differential capacitance and
diffusion in a half-space, respectively)

$\displaystyle H_{1}(\omega)$
$\displaystyle\sim-\frac{3}{\mathrm{i}\omega}-\frac{1}{5},$
$\displaystyle\text{as }\omega$ $\displaystyle\to 0,$ (47) $\displaystyle
H_{1}(\omega)$
$\displaystyle\sim\frac{1}{1-\sqrt{\mathrm{i}\omega}}\sim-\frac{1}{\sqrt{\mathrm{i}\omega}},$
$\displaystyle\text{as }\omega$ $\displaystyle\to\infty.$ (48)

Equation (43) for the fundamental is a commonly used impedance for batteries
(and many other electrochemical systems) that appears in linear EIS studies,
e.g. [21].

At second order, $O(\hat{I}^{2})$, there is a response at twice the excitation
frequency, $2\omega$, found to be

$\displaystyle z_{2,\pm}^{(2)}$
$\displaystyle=Z^{\mathrm{kin}}_{\mathrm{asym},\pm}+Z^{\mathrm{kin}}_{\mathrm{CD},\pm}+Z^{\mathrm{OCP}}_{\pm}+Z^{\mathrm{W}}_{2,\pm},$
(49)

with four contributions,

$\displaystyle Z^{\mathrm{kin}}_{\mathrm{asym},\pm}$
$\displaystyle=(\beta_{\pm}-1/2)\left(R_{\pm}^{(0)}\right)^{2},$ (50)
$\displaystyle Z^{\mathrm{kin}}_{\mathrm{CD},\pm}(\omega)$
$\displaystyle=R_{\pm}^{(0)\prime}\left(\frac{\xi_{\pm}}{D_{\pm}^{(0)}}\right)H_{1}\left(\frac{\omega}{D_{\pm}^{(0)}}\right),$
(51) $\displaystyle Z^{\mathrm{OCP}}_{\pm}(\omega)$
$\displaystyle=\frac{1}{2}U_{\pm}^{(0)\prime\prime}\left(\frac{\xi_{\pm}}{D_{\pm}^{(0)}}\right)^{2}\left[H_{1}\left(\frac{\omega}{D_{\pm}^{(0)}}\right)\right]^{2},$
(52) $\displaystyle Z^{\mathrm{W}}_{2,\pm}(\omega)$
$\displaystyle=U_{\pm}^{(0)\prime}\left(\frac{-D_{\pm}^{(0)\prime}}{D_{\pm}^{(0)}}\right)\left(\frac{\xi_{\pm}}{D_{\pm}^{(0)}}\right)^{2}H_{2}\left(\frac{\omega}{D_{\pm}^{(0)}}\right),$
(53)

These impedances, (50)-(53), each correspond to the different nonlinearities
in the model, listed in the Model Formulation section. The first and second
terms, (50)-(51), are due to nonlinearities in the kinetics. The kinetic term
$Z^{\mathrm{kin}}_{\mathrm{asym},\pm}$ represents the charge transfer
asymmetry, i.e. $\beta_{\pm}\neq 1/2$, and is purely real and independent of
frequency. The kinetic term $Z^{\mathrm{kin}}_{\mathrm{CD},\pm}$ is due to the
concentration dependence of the exchange current density (equivalently, the
charge transfer resistance $R_{\pm}(c_{\pm})$, see (24)), which is clear from
the presence of the derivative $R_{\pm}^{(0)\prime}$. The third term,
$Z^{\mathrm{OCP}}_{\pm}$, represents the nonlinearity of the OCP with respect
to stoichiometry, indicated by the second derivative appearing,
$U_{\pm}^{(0)\prime\prime}$. Finally, the term $Z^{\mathrm{W}}_{2,\pm}$ is due
to nonlinear diffusion of lithium within the particles, i.e. that the
diffusivity depends on concentration, indicated by $D_{\pm}^{(0)\prime}$.
Therefore, it could be interpreted as a second-order finite Warburg impedance,
cf. (44). Given our assumption on the functional dependence of $D_{\pm}^{(0)}$
(see (28)), we may substitute

$\displaystyle D_{\pm}^{(0)\prime}$
$\displaystyle=-\frac{1}{\tau_{\mathrm{d},\pm}}\left(U_{\pm}^{(0)\prime\prime}c_{\pm}^{(0)}+U_{\pm}^{(0)\prime}\right),$
(54)

into (53). At this order, the frequency dependence of the diffusion process is
encoded in $H_{2}(\omega)$, which is the solution of a linear ODE (see section
S2 of the supplementary material),found numerically. However, the high
frequency behaviour of $H_{2}$ may be extracted and is found to be

$\displaystyle H_{2}(\omega)$
$\displaystyle\sim\frac{1}{\mathrm{i}\omega}\left(1-\frac{1}{\sqrt{2}}\right),$
$\displaystyle\text{as }\omega$ $\displaystyle\to\infty,$ (55)

which is just a real multiple of $(H_{1})^{2}$ as $\omega\to\infty$ to leading
order. If the diffusion is linear, a common simplification, then
$D_{\pm}^{(0)\prime}=0$ and the term $Z^{\mathrm{W}}_{2,\pm}(\omega)$
vanishes.

In addition to the response at $2\omega$, there is a nonlinear response at
zero frequency, where many of the same terms arise:

$\displaystyle z_{0,\pm}^{(2)}$
$\displaystyle=2\mathrm{Re}\left[Z^{\mathrm{kin}}_{\mathrm{asym},\pm}+Z^{\mathrm{kin}}_{\mathrm{CD},\pm}+Z^{\mathrm{OCP}}_{0,\pm}+Z^{\mathrm{W}}_{0,\pm}\right],$
(56)

The terms are the same as in $z_{2,\pm}^{(2)}$, (49), except for the OCP and
Warburg terms, which are now

$\displaystyle Z^{\mathrm{OCP}}_{0,\pm}(\omega)$
$\displaystyle=\frac{1}{2}U_{\pm}^{(0)\prime\prime}\left(\frac{\xi_{\pm}}{D_{\pm}^{(0)}}\right)^{2}\left|H_{1}\left(\frac{\omega}{D_{\pm}^{(0)}}\right)\right|^{2},$
(57) $\displaystyle Z^{\mathrm{W}}_{0,\pm}(\omega)$
$\displaystyle=U_{\pm}^{(0)\prime}\left(\frac{-D_{\pm}^{(0)\prime}}{D_{\pm}^{(0)}}\right)\left(\frac{\xi_{\pm}}{D_{\pm}^{(0)}}\right)^{2}H_{0}\left(\frac{\omega}{D_{\pm}^{(0)}}\right),$
(58)

where a different function, $H_{0}(\omega)$ appears in the place of
$H_{2}(\omega)$—see section S2 of the supplementary material. Notice that
$z_{0,\pm}^{(2)}$ is real; this is because it represents a correction to the
time average of the response, and hence is not oscillatory.

The nonlinear impedances at $O(\hat{I}^{3})$ and higher can be readily
calculated, but with diminishing returns on their usability. Their amplitude
is inherently smaller, implying lower signal-to-noise ratios, with significant
increase in analytical complexity, so we do not consider them here.

The impedances for the full cell, including both electrodes and a series
resistance, are straightforwardly given in terms of the half-cell impedances
(43), (49), (56), using (22):

$\displaystyle z_{0}^{(0)}$ $\displaystyle=z_{0,+}^{(0)}-z_{0,-}^{(0)},$
$\displaystyle z_{1}^{(1)}$ $\displaystyle=z_{1,+}^{(1)}-z_{1,-}^{(1)}+R_{s},$
(59) $\displaystyle z_{2}^{(2)}$ $\displaystyle=z_{2,+}^{(2)}-z_{2,-}^{(2)},$
$\displaystyle z_{0}^{(2)}$ $\displaystyle=z_{0,+}^{(2)}-z_{0,-}^{(2)},$ (60)

Note that the series resistance $R_{s}$ only enters the $z_{1}^{(1)}$
expression as it is assumed to be a linear resistor.

The fundamental impedance $z_{1}^{(1)}$ is the usual linear one, measured
through conventional EIS, with dependence only on the series resistance
$R_{s}$ and electrode quantities $U_{+}^{(0)\prime}$, $U_{-}^{(0)\prime}$,
$D_{+}^{(0)}$, $D_{-}^{(0)}$, $\chi_{+}$, $\chi_{-}$, $\beta_{+}$,
$\beta_{-}$. The second order impedances, $z_{2}^{(2)}$ and $z_{0}^{(2)}$,
additionally depend on the higher derivatives $U_{+}^{(0)\prime\prime}$,
$U_{-}^{(0)\prime\prime}$, $D_{+}^{(0)\prime}$, $D_{-}^{(0)\prime}$. They also
have explicit dependence on $\beta_{+}$ and $\beta_{-}$ via the kinetic
asymmetry terms $Z^{\mathrm{kin}}_{\mathrm{asym},+}$ and
$Z^{\mathrm{kin}}_{\mathrm{asym},-}$.

###  3.3 SPM including double-layer capacitance

We now extend the formulae of the previous section to the more general case of
the SPM that includes double-layer capacitance, (15)-(24). One may repeat the
calculations of the previous section, starting from the full equations, but
the resulting impedance formulae may be expressed concisely in terms of those
already derived. The calculation of the half-cell impedances, now employing an
uppercase $Z$ to signify the inclusion of capacitance, is given in section
S1.2 of the supplementary material, and yields

$\displaystyle Z_{1,\pm}^{(1)}=\,\,$
$\displaystyle\frac{z_{1,\pm}^{(1)}(\omega)}{1\pm C_{\pm}\mathrm{i}\omega
z_{1,\pm}^{(1)}(\omega)}=\pm\left(\frac{1}{\left(C_{\pm}\mathrm{i}\omega\right)^{-1}}+\frac{1}{R_{\pm}^{(0)}+Z^{\mathrm{W}}_{1,\pm}(\omega)}\right)^{-1},$
(61) $\displaystyle Z_{2,\pm}^{(2)}=\,\,$
$\displaystyle\frac{z_{2,\pm}^{(2)}(\omega)}{\left[1\pm
C_{\pm}\mathrm{i}2\omega z_{1,\pm}^{(1)}(2\omega)\right]\left[1\pm
C_{\pm}\mathrm{i}\omega z_{1,\pm}^{(1)}(\omega)\right]^{2}},$ (62)
$\displaystyle Z_{0,\pm}^{(2)}=\,\,$
$\displaystyle\frac{z_{0,\pm}^{(2)}(\omega)}{\left|1\pm
C_{\pm}\mathrm{i}\omega z_{1,\pm}^{(1)}(\omega)\right|^{2}},$ (63)

where $z_{n,\pm}^{(m)}$ are the corresponding half-cell impedance formulae
assuming no capacitance ($C_{\pm}=0$), given by (43), (49) and (56). As there
is now more than one path the current can take through each electrode, these
formulae are no longer just a sum of individual impedances as in (43), (49)
and (56), but they are nonetheless still explicit.

The linear impedances $Z_{1,\pm}^{(1)}$ have the clear interpretation as being
the harmonic mean of $(C_{\pm}\mathrm{i}\omega)^{-1}$ and
$z_{1,\pm}^{(1)}(\omega)$. This corresponds to a capacitor in parallel with a
charge transfer resistance that is in series with a (finite) Warburg
element—this is exactly the Randles circuit. The nonlinear impedances do not
have such a clear analogy with traditional circuit elements. The expressions
(61)-(63) arise from nonlinearities in the charge transfer reaction and
diffusion process when they are placed in parallel with a capacitance, and
therefore had to be derived mathematically rather than from physical
intuition.

Figure 3: Equivalent circuit diagram representing the linear response of the
full cell model, consisting of a Randles circuit for each electrode along with
a series resistor $R_{s}$.

The impedances for a full-cell are given, using (22), by

$\displaystyle Z_{0}^{(0)}$ $\displaystyle=Z_{0,+}^{(0)}-Z_{0,-}^{(0)},$
$\displaystyle Z_{1}^{(1)}$ $\displaystyle=Z_{1,+}^{(1)}-Z_{1,-}^{(1)}+R_{s},$
(64) $\displaystyle Z_{2}^{(2)}$ $\displaystyle=Z_{2,+}^{(2)}-Z_{2,-}^{(2)},$
$\displaystyle Z_{0}^{(2)}$ $\displaystyle=Z_{0,+}^{(2)}-Z_{0,-}^{(2)}.$ (65)

The fundamental full-cell impedance then has the familiar equivalent circuit
analogy shown in Fig. 3, consisting of two Randles circuits (one for each
electrode) placed in series.

###  3.4 Simplifications: Separation of capacitance and diffusion timescales

Although the impedance formulae (61)-(63) are already explicit, some further
simplification is possible using common assumptions about the physical
timescales in the system. There are two prominent timescales in (61)-(63),
corresponding to capacitance and diffusion effects. These are given by
$R_{\pm}^{(0)}C_{\pm}$ and $1/D_{\pm}^{(0)}$, respectively. Note that these
timescales have been scaled by a known reference timescale $\tau^{*}$. Here
$1/D_{\pm}^{(0)}=\tau_{\mathrm{d},\pm}/(-U_{\pm}^{(0)\prime}c_{\pm}^{(0)})$
can be thought of as the _local_ diffusion timescale, at this DoD, which is
more convenient to use here instead of $\tau_{\mathrm{d},\pm}$ directly.
Typically the capacitance timescale ($O(10^{-2})$ s) is much shorter than the
diffusive one ($O(10^{4})$ s), hence we can consider the well-separated limit
of $R_{\pm}^{(0)}C_{\pm}\ll 1/D_{\pm}^{(0)}$.

Taking $R_{\pm}^{(0)}C_{\pm}\ll 1/D_{\pm}^{(0)}$ in (61)-(63) gives different
results depending on the frequency $\omega$. If the frequency is sufficiently
low, $\omega=O(D_{\pm}^{(0)})$ (region (I)), then diffusion will dominate and
the effects of capacitance will be negligible, i.e. $C_{\pm}$ may be set to
zero. If $\omega=O(D_{\pm}^{(0)}/(R_{\pm}^{(0)}C_{\pm}))$ (region (II)), the
opposite is true and we may take the high $\omega$ limits of
$z_{n,\pm}^{(m)}(\omega)$, which are real and given by

$\displaystyle z_{1,\pm}^{(1)}(\omega)$ $\displaystyle\to\pm
R_{\pm}^{(0)},\qquad\omega\to\infty,$ (66) $\displaystyle
z_{2,\pm}^{(2)}(\omega)$
$\displaystyle\to(\beta_{\pm}-1/2)(R_{\pm}^{(0)})^{2},\qquad\omega\to\infty,$
(67) $\displaystyle z_{0,\pm}^{(2)}(\omega)$ $\displaystyle\to
2(\beta_{\pm}-1/2)(R_{\pm}^{(0)})^{2},\qquad\omega\to\infty.$ (68)

Then approximations of (61)-(63) in both regions take the form

$\displaystyle Z_{1,\pm}^{(1)}\sim\,\,$
$\displaystyle\begin{dcases}z_{1,\pm}^{(1)}(\omega)&\text{(I),
}\omega=O(D_{\pm}^{(0)}),\\\ \frac{\pm
R_{\pm}^{(0)}}{1+R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega}&\text{(II),
}\omega=O\left(\frac{D_{\pm}^{(0)}}{R_{\pm}^{(0)}C_{\pm}}\right),\\\
\end{dcases}$ (69) $\displaystyle Z_{2,\pm}^{(2)}\sim\,\,$
$\displaystyle\begin{dcases}z_{2,\pm}^{(2)}(\omega)&\text{(I),
}\omega=O(D_{\pm}^{(0)}),\\\
\frac{(\beta_{\pm}-1/2)(R_{\pm}^{(0)})^{2}}{\left[1+2R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega\right]\left[1+R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega\right]^{2}},&\text{(II),
}\omega=O\left(\frac{D_{\pm}^{(0)}}{R_{\pm}^{(0)}C_{\pm}}\right),\\\
\end{dcases}$ (70) $\displaystyle Z_{0,\pm}^{(2)}\sim\,\,$
$\displaystyle\begin{dcases}z_{0,\pm}^{(2)}(\omega)&\text{(I),
}\omega=O(D_{\pm}^{(0)}),\\\
\frac{2(\beta_{\pm}-1/2)(R_{\pm}^{(0)})^{2}}{\left|1+R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega\right|^{2}},&\text{(II),
}\omega=O\left(\frac{D_{\pm}^{(0)}}{R_{\pm}^{(0)}C_{\pm}}\right),\\\
\end{dcases}$ (71)

Figure 4: Equivalent circuit diagrams representing the linear response of each
electrode, but corresponding to different impedance formulae. (a) Formula
(61), with no assumption on the diffusion and capacitance timescales. (b)
Composite (approximate) formula (73), where the timescales are assumed well-
separated and the Warburg impedance acts in series with the kinetics.

There is a narrow region between Region I and region II where both capacitance
and diffusive effects act with similar small strength. This narrow region
occurs for frequencies

$\displaystyle\omega$
$\displaystyle\approx\left(\frac{-{U_{\pm}^{(0)}}^{\prime}\xi_{\pm}\sqrt{D_{\pm}^{(0)}}}{\left(R_{\pm}^{(0)}\right)^{2}C_{\pm}}\right)^{2/3}.$
(72)

We can now create an approximation valid for all the frequencies by forming a
composite approximation made by summing the region I and II approximations and
subtracting the leading order solution in the overlap region. This overlap
solution is found by taking the $\omega\to\infty$ limit from region I (see
(66)-(68)) or the $\omega\to 0$ limit from region II (at the fundamental
frequency, this is simply the charge transfer resistance $R_{\pm}^{(0)}$). The
resulting composite solutions are

$\displaystyle Z_{1,\pm}^{(1),\mathrm{comp}}$
$\displaystyle\sim\pm\left[\frac{R_{\pm}^{(0)}}{1+R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega}+Z^{\mathrm{W}}_{1,\pm}(\omega)\right],$
(73) $\displaystyle Z_{2,\pm}^{(2),\mathrm{comp}}$
$\displaystyle\sim\frac{Z^{\mathrm{kin}}_{\mathrm{asym},\pm}}{\left[1+2R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega\right]\left[1+R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega\right]^{2}}+Z^{\mathrm{kin}}_{\mathrm{CD},\pm}+Z^{\mathrm{OCP}}_{\pm}+Z^{\mathrm{W}}_{2,\pm},$
(74) $\displaystyle Z_{0,\pm}^{(2),\mathrm{comp}}$
$\displaystyle\sim\frac{2Z^{\mathrm{kin}}_{\mathrm{asym},\pm}}{\left|1+R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega\right|^{2}}+2\mathrm{Re}\left[Z^{\mathrm{kin}}_{\mathrm{CD},\pm}+Z^{\mathrm{OCP}}_{0,\pm}+Z^{\mathrm{W}}_{0,\pm}\right].$
(75)

Notice that these are close in form to the expressions (43), (49) and (56),
but with capacitance modifying only the first term. Looking at the linear
impedance (73), it now appears that the Warburg term is no longer in
_parallel_ with the capacitance (as in Fig. 4$(a)$) but in _series_ (as in
Fig. 4$(b)$). Since they each dominate at different frequency scales, the
diffusion process interacts minimally with the kinetics. The same is true for
the nonlinear impedances (74)-(75), where the terms that depend on the
diffusion process separate out, leaving only the constant kinetic term
$Z^{\mathrm{kin}}_{\mathrm{asym},\pm}$ (due to reaction asymmetry) in
“parallel” with the capacitance. This simplifies the structure and
interpretation of impedances, since the contribution of each nonlinearity can
be more easily assessed.

##  4 Methods

Having introduced the model and analysed it in the frequency domain, we now
detail the numerical, experimental and parameterization methods employed in
the remainder of the paper to produce the results given in the Results
section.

###  4.1 Numerical solution in the time domain for generation of synthetic
data

Here we describe the numerical methods for solving the SPM in the time domain
and the procedure for generating synthetic impedance data from a known set of
parameter values—this is useful for parameter identifiability and error
sensitivity analysis (the Results section) because ‘ground truth’ parameters
are known exactly. Synthetic data was generated by solving the full nonlinear
model equations (15)-(24) in the time domain numerically given a sinusoidal
input current of dimensional frequency $\omega^{*}$ and amplitude
$2\hat{I}^{*}$ (or peak-to-peak amplitude $4\hat{I}^{*}$). This was then
converted to the frequency domain with a Fourier transform. The radial
dimension within the electrode particles was discretized using a finite volume
scheme, and the time integration was performed using the adaptive explicit
stiff ODE solver ode15s in MATLAB. To ensure sufficient relaxation of initial
transients and convergence of the output voltage to a periodic profile, 20
periods were simulated and only the final 2 were selected to be fast Fourier
transformed (FFT) with 20 equispaced samples per period. This procedure was
performed for 30 frequencies logarithmically spaced from $\omega^{*}=$
${10}^{-4}\text{\,}\mathrm{H}\mathrm{z}$ to
${10}^{2}\text{\,}\mathrm{H}\mathrm{z}$, i.e. 5 frequencies per decade. The
parameter set used corresponds to a LCO/LiC6 cell, with values modified from
[26, 15], summarized in Tables S1, S2 of the supplementary material.
Simulations were performed at three DoDs (30%, 50% and 70%) with a fixed
current amplitude of $50\text{\,}\mathrm{m}\mathrm{A}$.

Zero-mean Gaussian noise was added to the simulated voltages in the frequency
domain ($n=1,2$)

$\displaystyle\widehat{V}_{n}^{*,\mathrm{data}}$
$\displaystyle=\widehat{V}_{n}^{*}+\varepsilon_{n}^{*},$
$\displaystyle\varepsilon_{n}^{*}$
$\displaystyle=\varepsilon_{n,r}^{*}+\mathrm{i}\varepsilon_{n,i}^{*},\quad\mbox{where
}\varepsilon_{n,r}^{*},\varepsilon_{n,i}^{*}\sim N(0,\sigma_{n}^{*2}),$ (76)

and the standard deviations were chosen to give similar signal-to-noise ratios
to those seen in our experiments, i.e. $\sigma_{1}^{*}=\sigma_{2}^{*}=10^{-7}$
V. The synthetic impedances were calculated from

$\displaystyle Z_{n}^{*,\mathrm{data}}$
$\displaystyle=\frac{\widehat{V}_{n}^{*,\mathrm{data}}}{(\hat{I}^{*})^{n}},\qquad
n=1,2.$ (77)

###  4.2 Experimental measurements

Experiments were undertaken using a bespoke test rig to create an NLEIS
dataset to parameterize the model. All measurements were collected in the time
domain, then post-processed in MATLAB. A bipolar power source/sink (Kikusui
PBZ60-6.7) was connected in series with an NMC/graphite pouch cell (Kokam
SLPB533459H4, $740\text{\,}\mathrm{m}\mathrm{A}\mathrm{h}$ [27]) and a
$0.1\text{\,}\mathrm{\SIUnitSymbolOhm}$ precision shunt resistor as shown in
Fig. S1 in the supplementary material. Voltage measurements were made across
the test battery and the shunt respectively using a 24-bit analog-to-digital
converter (National Instruments 9239). Two cells were tested.

Prior to NLEIS experiments, the discharge capacity of each cell was measured
at $19\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and found to be
$720\text{\,}\mathrm{m}\mathrm{A}\mathrm{h}$ and
$729\text{\,}\mathrm{m}\mathrm{A}\mathrm{h}$ respectively. For capacity
measurements, first a constant-current charge of
$500\text{\,}\mathrm{m}\mathrm{A}$ to $4.2\text{\,}\mathrm{V}$ was applied
followed by a constant-voltage charge until the current dropped below
$50\text{\,}\mathrm{m}\mathrm{A}$. Then, after a $30\text{\,}\mathrm{min}$
rest period, a constant-current discharge at
$500\text{\,}\mathrm{m}\mathrm{A}$ to a voltage of $2.7\text{\,}\mathrm{V}$
was conducted followed by a constant-voltage discharge to a current below
$50\text{\,}\mathrm{m}\mathrm{A}$.

Nonlinear EIS measurements were collected at 10% DoD to 90% DoD (equivalently
90% SOC to 10% SOC), in steps of 10%. At each DoD, zero-mean sinusoidal
currents of $1\text{\,}\mathrm{A}$ peak-to-peak were applied and voltage
measurements taken at frequencies from $0.1\text{\,}\mathrm{H}\mathrm{z}$ to
$10\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$, with 10 measurements per decade.
The length of the measurement was taken as the larger of
$2\text{\,}\mathrm{s}$ or 20 cycles at each frequency. For more details on the
procedures, including thermal limits, data acquisition and processing, see
section S4 of the supplementary material.

The measured voltage and current signals were each multiplied by a Hann window
function, and then fast-Fourier transformed using MATLAB’s built-in FFT
function, giving the components $\widehat{V}_{n}(\omega,\hat{I})$ of the
voltage as per (35).

###  4.3 OCP functions and electrode balancing

Individual electrode OCPs $U_{\pm}^{*}$ for the Kokam SLPB533459H4 cells were
determined from 3-electrode measurements by McTurk et al. [25], who inserted a
minimally invasive lithium reference electrode into cells of the same model.
This data for $U_{\pm}^{*}$ is given in terms of normalised discharge capacity
$Q=\text{DoD}/100$ with $Q\in[0,1]$. For use with our model, analytical
expressions (S4.5 of supplementary material) were fitted to the data, shown in
Figs. 5($a$)-($b$). Fig. 5($c$) demonstrates that the resulting full-cell OCV
agrees with the measured OCV of the cells used in this paper.

Figure 5: $(a)$ Negative electrode OCP, $(b)$ positive electrode OCP, and
$(c)$ full-cell OCV for a Kokam SLPB533459H4
$740\text{\,}\mathrm{m}\mathrm{A}\mathrm{h}$ pouch cell. Dot markers are
3-electrode data from [25], solid lines are analytical expressions (Eqs.
(S4.1), (S4.2) in supplementary material) that were fitted to individual
electrode OCP data. Panel $(c)$ compares these fits to the measured full-cell
OCV for the present cells from which NLEIS data was collected.

To balance the electrodes, the stoichiometry limits $c_{\pm}^{0\%}$ and
$c_{\pm}^{100\%}$ corresponding to 0% and 100% DoD ($Q=0$ and 1) were chosen
using knowledge of the electrode chemistries. The resulting stoichiometry
limits (and values of $\xi_{\pm}$ that follow from (33)) are given later in
Table 1. Then $U_{\pm}^{*}$ is expressed in terms of $c_{\pm}$ using (32).
When it is required, the first derivative of $U_{\pm}^{*}$ is calculated from
the analytical fits, but the second derivative is calculated numerically from
the OCP data. The OCP expressions and further details of the fitting
procedures and electrode balancing are given in section S4.5 of the
supplementary material.

###  4.4 Parameter estimation

To identify parameters and analyse their uniqueness, we consider fitting the
(dimensional) model impedances $Z_{n}^{*,(n)}=\Phi^{*}Z_{n}^{(n)}/(J^{*})^{n}$
to measured or simulated data $Z_{n}^{*,\mathrm{data}}$ across a range of
frequencies at a given DoD. The parameter estimation problem takes the form

$\displaystyle\hat{\bm{\theta}}$
$\displaystyle=\underset{\bm{\theta}}{\mathrm{argmin}}\,l(\bm{\theta}),$ (78)

where $\bm{\theta}$ is the vector of parameters to fit, and $l$ is a loss
function describing the goodness of fit. It is not obvious what loss function
to use for NLEIS data since more than one impedance (i.e. $Z_{1}^{*}$ and
$Z_{2}^{*}$) is fitted simultaneously, and higher harmonics are typically
smaller in magnitude but also more sensitive to some parameters. A
straightforward choice of the total squared error across both harmonics would
not account for this difference in scale or magnitude. Therefore, the relative
weighting of the error across each harmonic merits consideration.

We take the approach of maximum likelihood estimation (MLE), assuming zero-
mean Gaussian measurement noise added to the voltage harmonics, as in (76).
Notably, we assume in general that the noise variances of the first and second
harmonics are not the same. The details can be found in section S5 of the
supplementary material, and the resulting MLE for $\bm{\theta}$ given NLEIS
data $Z_{1,j}^{*,\mathrm{data}}$, $Z_{2,j}^{*,\mathrm{data}}$ at frequencies
$\omega_{j},j=1,2,\ldots,N_{\omega}$, is reduced to

$\displaystyle\hat{\bm{\theta}}_{\mathrm{NLEIS}}$
$\displaystyle=\underset{\bm{\theta}}{\mathrm{argmin}}\,(l_{1}(\bm{\theta})+l_{2}(\bm{\theta})),$
(79)

where

$\displaystyle l_{1}(\bm{\theta})$
$\displaystyle=\log\left(\sum_{j=1}^{N_{\omega}}\left|Z_{1}^{*(1)}(\omega_{j};\bm{\theta})-Z_{1,j}^{*,\mathrm{data}}\right|^{2}\right),$
(80) $\displaystyle l_{2}(\bm{\theta})$
$\displaystyle=\log\left(\sum_{j=1}^{N_{\omega}}\left|Z_{2}^{*(2)}(\omega_{j};\bm{\theta})-Z_{2,j}^{*,\mathrm{data}}\right|^{2}\right),$
(81)

are the (negative) log-likelihoods, up to a scaling factor, from each
harmonic. We will refer to their sum as the total loglikelihood
$l_{12}(\bm{\theta})=l_{1}+l_{2}$. The form of $l_{12}$, which is the sum of
the logarithm of the squared error of each harmonic, originates from the MLE
process and has the advantage of being independent of the scale of each
harmonic, as desired. That is, rescaling
$(Z_{n}^{*n},Z_{n}^{*,\mathrm{data}})\mapsto(aZ_{n}^{*n},aZ_{n}^{*,\mathrm{data}})$
for any $a$ only changes $l_{12}$ by an additive constant, leaving
$\mathrm{argmin}_{\bm{\theta}}\>l_{12}$ unchanged. Given only linear EIS data,
consisting of $Z_{1,j}^{*,\mathrm{data}}$, the corresponding MLE is

$\displaystyle\hat{\bm{\theta}}_{\mathrm{EIS}}$
$\displaystyle=\underset{\bm{\theta}}{\mathrm{argmin}}\,l_{1}(\bm{\theta}),$
(82)

where the loss function is simply $l_{1}$.

The minimization problem (78) was solved using nonlinear optimization routines
(local and global) in MATLAB employing bound constraints as described in
section S5 of the supplementary material.

##  5 Results

###  5.1 Analysis of model impedances

We begin with an illustration of the form of the exact model impedances, up to
second order, given by the formulae (61), (62) and (63). The dependence of
$Z_{2,\pm}^{(2)}$ in particular on the various model parameters will be
demonstrated. Then, the simplified (composite) formulae (73), (74) and (75)
will be compared to the aforementioned exact results.

Typical Nyquist plots of the impedances $Z_{1,+}^{(1)}$, $Z_{2,+}^{(2)}$ and
$Z_{0,+}^{(2)}$ are shown in Fig. 6, using the parameter values in Tables S1,
S2 in the supplementary material. (Only the positive electrode impedances are
shown—negative electrode ones have a similar structure.) In Fig. 6$(a)$,
$Z_{1,+}^{(1)}$ consists of a high frequency semi-circular kinetic arc, and a
low frequency diffusive tail with a capacitive effect
($Z_{1,+}^{(1)}=O(1/\omega)$) as $\omega\to 0$. That is, the same as the well-
understood Randles circuit with finite-space diffusion. In Fig. 6$(b)$, we see
the second harmonic impedance $Z_{2,+}^{(2)}$ also consists of a high
frequency kinetic arc and diffusive tail, but this arc appears “spiral” in
nature rather than circular. Equation (70) (region II) describes this and,
excluding the factor of $(1+2R_{\pm}^{(0)}C_{\pm}\mathrm{i}\omega)$ (which
affects mainly the shape close to the origin), shows that $Z_{2,+}^{(2)}$ is
proportional to the square of the semi-circular kinetics in $Z_{1,\pm}^{(1)}$,
resulting in a _cardioid_ -like shape. Lastly, the second-order correction at
the zeroth harmonic, $Z_{0,+}^{(2)}$, contains similar model information as
$Z_{2,+}^{(2)}$ but with no imaginary component. In Fig. 6$(c)$, we see
$Z_{0,+}^{(2)}=O(\omega^{-2})$ at high frequencies (square of double-layer
capacitive effects) but also begin to see this behaviour at low frequencies
(square of electrode “differential capacitive” effects). However, in
intermediate regions $Z_{0,+}^{(2)}$ appears relatively flat and the various
effects are difficult to delineate.

Next, we demonstrate in Fig. 7 how the various features of $Z_{2,+}^{(2)}$
(and therefore $Z_{2,-}^{(2)}$) depend on quantities in the model. Starting
from the values in Tables S1, S2 and at 30% DoD, we varied $R_{+}^{(0)}$,
$C_{+}$, $R_{+}^{(0)\prime}$, $\beta_{+}$, $\tau_{\mathrm{d},+}$, and
$U_{+}^{(0)\prime\prime}$ across a representative range (considering both
positive and negative values for the derivatives $R_{+}^{(0)\prime}$,
$U_{+}^{(0)\prime\prime}$). From Fig. 7 we observe that:

  1. $(a)$

$R_{+}^{(0)}$ affects the size of the kinetic arc (its “width” on the real
axis is
$Z^{\mathrm{kin}}_{\mathrm{asym},+}=(\beta_{+}-1/2)(R_{+}^{(0)})^{2}$);

  2. $(b)$

$C_{+}$ affects the frequency dependence of the kinetic arc, but not its shape
(just as for $Z_{1,+}^{(1)}$);

  3. $(c)$

$R_{+}^{(0)\prime}$ affects the initial direction and magnitude of the
diffusion tail (due to $Z^{\mathrm{kin}}_{\mathrm{CD},+}$) emanating from the
kinetic arc;

  4. $(d)$

$\beta_{+}$ affects the size (via the magnitude $|\beta_{+}-1/2|$) and
orientation (via $\mathrm{sgn}(\beta_{+}-1/2$)) of the kinetic arc (see point
$(a)$ above). Important to note that as $\beta_{+}\to 1/2$ the totality of the
arc shrinks to the origin. Thus $\beta_{+}\neq 1/2$ is _necessary_ to observe
an arc at all (this was first pointed out by Murbach et al. [5]);

  5. $(e)$

$\tau_{\mathrm{d},+}$ affects the length and shape of the diffusion tail (via
the terms $Z^{\mathrm{kin}}_{\mathrm{CD},+}$ and $Z^{\mathrm{OCP}}_{+}$);

  6. $(f)$

$U_{+}^{(0)\prime\prime}$ affects whether the diffusion tail ultimately
diverges to $\mathrm{Re}=-\infty$ (when $U_{+}^{(0)\prime\prime}>0$) or
$\mathrm{Re}=+\infty$ (when $U_{+}^{(0)\prime\prime}<0$) as $\omega\to 0$.

It is clear that the possible behaviour of $Z_{2,+}^{(2)}$ is much more varied
than $Z_{1,+}^{(1)}$, with visible signatures of many quantities in the model
that are not visible with EIS, e.g. $\beta_{+}$, $R_{+}^{(0)\prime}$, and
$U_{+}^{(0)\prime\prime}$. In addition, dependence of $Z_{2,+}^{(2)}$ on
$R_{+}^{(0)}$, $C_{+}$ and $\tau_{\mathrm{d},+}$ gives further information
that may be used to improve their identifiability over traditional linear EIS.

Figure 6: Nyquist plot of the fundamental (linear) impedance $Z_{1,+}^{(1)}$,
and the second-order nonlinear impedances $Z_{2,+}^{(2)}$ (second harmonic)
and $Z_{0,+}^{(2)}$ (zero frequency mode), as calculated from formulae (61),
(62) and (63). Parameters in Table S2 at $\mathrm{DoD}=30\%$ (i.e.
$c_{+}=0.80$). Frequency range is $10^{-4}\leq\omega/2\pi\leq 10^{2}$ [Hz],
with 5 values (marked) per decade. Only the impedances for positive electrode
shown.
Figure 7: Dependence of the second harmonic impedance $Z_{2,+}^{(2)}$ on
quantities $(a)$ $R_{+}^{(0)}$, $(b)$ $C_{+}$, $(c)$ $R_{+}^{(0)\prime}$,
$(d)$ $\beta_{+}$, $(e)$ $\tau_{\mathrm{d},+}$, $(f)$
$U_{+}^{(0)\prime\prime}$. For each, we used the value in Table S2 (at
$\mathrm{DoD}=30\%$), shown in blue, and 2 other representative values.

Finally, we remark on the accuracy of the simplified (composite) expressions
(73), (74) and (75), which assume that the capacitive timescale
$R_{\pm}^{(0)}C_{\pm}$ is much shorter than the diffusive one,
$\tau_{\mathrm{d},\pm}=1/D_{\pm}^{(0)}$. Comparing these to the exact
impedances for a full-cell, we find them practically indistinguishable—across
all frequencies and each harmonic—with mean relative errors less than 0.03%
for the fundamental and 0.64% for the second harmonic. Details are given in
section S6 of the supplementary material. Consequently, these composite
expressions are a useful and simple substitute for the full expressions in
practice.

###  5.2 Structural parameter identifiability

Nine dynamical parameter groups (34) are sufficient to fully parameterize the
nonlinear SPM (15)-(24). However, these parameters are not necessarily able to
be determined uniquely (or at all) from the model structure. This problem is
referred to as structural parameter identifiability [28], and we will consider
the improvements that NLEIS provides relative to linear EIS [21].

Structural identifiability is typically defined for linear systems, but here
we will use a nonlinear extension. Consider a set of nonlinear impedances
$\\{Z_{1}^{(1)}(\omega,\bm{\theta}),Z_{2}^{(2)}(\omega,\bm{\theta}),\ldots\\}$,
characterising a nonlinear model in the frequency domain, for parameters
$\bm{\theta}\in\Theta$. If the system of equations with
$\bm{\theta},\tilde{\bm{\theta}}\in\Theta$

$\displaystyle Z_{1}^{(1)}(\omega,\bm{\theta})$
$\displaystyle=Z_{1}^{(1)}(\omega,\tilde{\bm{\theta}}),$ (83) $\displaystyle
Z_{2}^{(2)}(\omega,\bm{\theta})$
$\displaystyle=Z_{2}^{(2)}(\omega,\tilde{\bm{\theta}}),$ (84)
$\displaystyle\vdots$

for all $\omega\in\mathbb{R}$, has: (i) a unique solution for $\bm{\theta}$,
then the model is globally identifiable; (ii) a finite number of solutions for
$\bm{\theta}$, it is locally identifiable; (iii) an infinite number of
solutions for $\bm{\theta}$, it is unidentifiable.

####  5.2.1 Single electrode

First, we will consider the parameter identifiability of a single electrode
(positive or negative) and restrict ourselves to the accurate simplified
formulae (73)-(74). For a single ($\pm$) electrode, the relevant unknown
parameters are $\chi_{\pm},C_{\pm},\beta_{\pm},\tau_{\mathrm{d},\pm}$, but for
convenience we will consider

$\displaystyle R_{\pm}^{(0)}$
$\displaystyle=\frac{2\chi_{\pm}}{(c_{\pm}^{(0)})^{\beta_{\pm}}(1-c_{\pm}^{(0)})^{1-\beta_{\pm}}},$
(85)

in place of $\chi_{\pm}$, and discuss identifiability of $\chi_{\pm}$
subsequently.

For a single electrode, if only $Z_{1,\pm}^{(1)}$ is given, (73) shows that
$R_{\pm}^{(0)},C_{\pm}$ have unique solutions (globally identifiable), and so
does $\tau_{\mathrm{d},\pm}$ provided $U_{\pm}^{(0)\prime}\neq 0$. If
$U_{\pm}^{(0)\prime}=0$, i.e. if the OCP slope is exactly zero, the diffusive
contribution to the impedance vanishes and $\tau_{\mathrm{d},\pm}$ is
unidentifiable—as per Bizeray et al. [21] for Fickian diffusion. However, the
charge transfer coefficient $\beta_{\pm}$ does not appear in $Z_{1,\pm}^{(1)}$
and is therefore always unidentifiable. In this case $\chi_{\pm}$ cannot be
identified from (85), leaving only $R_{\pm}^{(0)}$ identified.

However, if the second harmonic is also considered, i.e.,
$\\{Z_{1,\pm}^{(1)},Z_{2,\pm}^{(2)}\\}$ is given, then more information on
each parameter is provided. In particular, $\beta_{\pm}$ (and hence
$\chi_{\pm}$) becomes globally identifiable since it appears explicitly in
$Z_{2,\pm}^{(2)}$. Regarding the diffusion timescale $\tau_{\mathrm{d},\pm}$:
if $U_{\pm}^{(0)\prime}=0$ (hence $\tau_{\mathrm{d},\pm}$ unidentifiable from
$Z_{1,\pm}^{(1)}$ alone) but $U_{\pm}^{(0)\prime\prime}\neq 0$, then the terms
$Z_{\pm}^{\mathrm{OCP}}$ and $Z_{\pm,2}^{\mathrm{W}}$ in $Z_{2,\pm}^{(2)}$ do
not vanish and in fact make $\tau_{\mathrm{d},\pm}$ identifiable. However, if
$U_{\pm}^{(0)\prime\prime}=0$ also, then $\tau_{\mathrm{d},\pm}$ remains
unidentifiable.

####  5.2.2 Full two-electrode cell

Given the impedances $\\{Z_{1,\pm}^{(1)},Z_{2,\pm}^{(2)}\\}$ for a single
electrode, the full cell impedances are given by (64)-(65) and they inherit
each electrode’s structural identifiability since no parameters appear in
impedances from _both_ electrodes. However, they may have additional
identifiability issues arising from cell symmetries, i.e., “electrode
swapping”, or inability to decide whether a parameter refers to the positive
or negative electrode.

If the full-cell linear impedance
$Z_{1}^{(1)}=Z_{1,+}^{(1)}-Z_{1,-}^{(1)}+R_{s}$ is given, the series
resistance $R_{s}$ is unique (globally identifiable) and can be extracted from
the high frequency limit $\omega\to\infty$, where $Z_{1}^{(1)}\to R_{s}$.
However, the identifiable kinetic parameters now have two solutions: if
$(R_{-}^{(0)},C_{-},R_{+}^{(0)},C_{+})$ is a solution then so is
$(R_{+}^{(0)},C_{+},R_{-}^{(0)},C_{-})$. Hence they are only _locally_
identifiable. A similar symmetry exists for the diffusion timescales
$\tau_{\mathrm{d},\pm}$ but only for special values of the OCP slopes.
Substituting (73) into condition (83), and keeping only the Warburg terms,
gives (recall
$D_{\pm}^{(0)}=-U_{\pm}^{(0)\prime}c_{\pm}^{(0)}/\tau_{\mathrm{d,\pm}}$):

$\displaystyle\frac{\xi_{+}U_{+}^{(0)\prime}}{D_{+}^{(0)}}H_{1}\left(\frac{\omega}{D_{+}^{(0)}}\right)+\frac{\xi_{-}U_{-}^{(0)\prime}}{D_{-}^{(0)}}H_{1}\left(\frac{\omega}{D_{-}^{(0)}}\right)$
$\displaystyle=\frac{\xi_{+}U_{+}^{(0)\prime}}{\tilde{D}_{+}^{(0)}}H_{1}\left(\frac{\omega}{\tilde{D}_{+}^{(0)}}\right)+\frac{\xi_{-}U_{-}^{(0)\prime}}{\tilde{D}_{-}^{(0)}}H_{1}\left(\frac{\omega}{\tilde{D}_{-}^{(0)}}\right).$
(86)

If the OCP slopes are nonzero and satisfy
$\xi_{+}U_{+}^{(0)\prime}=\xi_{-}U_{-}^{(0)\prime}$, this reduces to

$\displaystyle\frac{1}{D_{+}^{(0)}}H_{1}\left(\frac{\omega}{D_{+}^{(0)}}\right)+\frac{1}{D_{-}^{(0)}}H_{1}\left(\frac{\omega}{D_{-}^{(0)}}\right)$
$\displaystyle=\frac{1}{\tilde{D}_{+}^{(0)}}H_{1}\left(\frac{\omega}{\tilde{D}_{+}^{(0)}}\right)+\frac{1}{\tilde{D}_{-}^{(0)}}H_{1}\left(\frac{\omega}{\tilde{D}_{-}^{(0)}}\right),$
(87)

The left hand side has the symmetry $D_{-}^{(0)}\mapsto D_{+}^{(0)}$,
$D_{+}^{(0)}\mapsto D_{-}^{(0)}$, meaning there are always two solutions for
$(D_{+}^{(0)},D_{-}^{(0)})$ and hence
$(\tau_{\mathrm{d,+}},\tau_{\mathrm{d,-}})$. This symmetry in terms of
$\tau_{\mathrm{d,\pm}}$ is

$\displaystyle\tau_{\mathrm{d,+}}$
$\displaystyle\mapsto\frac{U_{+}^{(0)\prime}c_{+}^{(0)}}{U_{-}^{(0)\prime}c_{-}^{(0)}}\tau_{\mathrm{d,-}},$
$\displaystyle\tau_{\mathrm{d,-}}$
$\displaystyle\mapsto\frac{U_{-}^{(0)\prime}c_{-}^{(0)}}{U_{+}^{(0)\prime}c_{+}^{(0)}}\tau_{\mathrm{d,+}}.$
(88)

This electrode symmetry of diffusion timescales in $Z_{1}^{(1)}$ was pointed
out by Bizeray _et al._ [21], which we extend here to the case of nonlinear
diffusion.

If the full-cell second harmonic is also considered, i.e.,
$\\{Z_{1,}^{(1)},Z_{2}^{(2)}\\}$ is given, the above symmetries are removed in
most cases. Substituting (74) into condition (84), as $H_{1}$, $H_{1}^{2}$ and
$H_{2}$ are distinct functions of $\omega$, the diffusion symmetry
$D_{-}^{(0)}\mapsto D_{+}^{(0)}$, $D_{+}^{(0)}\mapsto D_{-}^{(0)}$ only holds
(using similar arguments to (86)-(87)) if
$\xi_{+}^{2}U_{+}^{(0)\prime\prime}=-\xi_{-}^{2}U_{-}^{(0)\prime\prime}$ and
$\xi_{+}^{2}U_{+}^{(0)\prime}/c_{+}^{(0)}=-\xi_{-}^{2}U_{-}^{(0)\prime}/c_{-}^{(0)}$.
But we already require $\xi_{+}U_{+}^{(0)\prime}=\xi_{-}U_{-}^{(0)\prime}$
from $Z_{1}^{(1)}$ which, if $U_{\pm}^{(0)\prime}\not=0$, reduces the latter
to $\xi_{+}/c_{+}^{(0)}=-\xi_{-}/c_{-}^{(0)}$, and this is not possible since
$\xi_{+}$ and $\xi_{-}$ are both positive. If actually
$U_{\pm}^{(0)\prime}=0$, then the symmetry can persist so long as
$\xi_{+}^{2}U_{+}^{(0)\prime\prime}=-\xi_{-}^{2}U_{-}^{(0)\prime\prime}$.

If the diffusion symmetry is removed, then the terms
$Z_{\pm,\mathrm{CD}}^{\mathrm{kin}}$ (see (51)) prevent the symmetry in
$R_{\pm}^{(0)},C_{\pm}$ too as they depend on both kinetic and diffusive
parameters. However, in the special case of
$R_{-}^{(0)\prime}=R_{+}^{(0)\prime}=0$, i.e. both exchange current densities
are independent of concentration (at this DoD), then
$Z_{\pm,\mathrm{CD}}^{\mathrm{kin}}$ vanishes and a symmetry in the kinetic
parameters does persist, given by $R_{\pm}^{(0)}\mapsto R_{\mp}^{(0)}$,
$C_{\pm}\mapsto C_{\mp},\beta_{\pm}\mapsto 1-\beta_{\mp}$.

In summary, including the second harmonic $Z_{2}^{(2)}$, makes $\beta_{\pm}$
(and hence $\chi_{\pm}$) globally identifiable, and improves the
identifiability of the diffusion timescales $\tau_{\mathrm{d},\pm}$ when an
OCP is flat ($U_{\pm}^{(0)\prime}=0$) but has nonzero curvature
($U_{\pm}^{(0)\prime\prime}\neq 0$). In addition, it removes the electrode
swapping symmetries in general, except when both OCVs are flat
($U_{+}^{(0)\prime}=U_{-}^{(0)\prime}=0$) or both exchange current densities
are independent of concentration ($R_{-}^{(0)\prime}=R_{+}^{(0)\prime}=0$).

###  5.3 Parameter estimation from synthetic data

With the structural identifiability from the first two harmonics analysed, we
move to analyse the practical identifiability, i.e. parameter estimation from
data. In this section, we consider noisy synthetic data from a known parameter
set—see the Methods section for the generation method and parameter values.
The parameter estimation algorithm from EIS data ($Z_{1}^{*,\mathrm{data}}$
only) and NLEIS data ($Z_{1}^{*,\mathrm{data}}$ and $Z_{2}^{*,\mathrm{data}}$)
was given in the Methods section and results in estimators
$\hat{\bm{\theta}}_{\mathrm{EIS}}$, $\hat{\bm{\theta}}_{\mathrm{NLEIS}}$ for
the parameter groups.

To test the consistency and noise sensitivity of our approach we repeatedly
applied the estimation algorithm to synthetic data but each time with separate
samples drawn from the noise distribution (given in (76)). This resampling,
known as “bootstrapping” [29], allows the construction of the distribution of
the estimators $\hat{\bm{\theta}}_{\mathrm{EIS}}$,
$\hat{\bm{\theta}}_{\mathrm{NLEIS}}$. The distributions (histograms) and their
means for 200 resamples are shown and compared in Figs. 8 (with data at 50%
DoD) and 9 (at 70% DoD). Note that the transfer coefficient $\beta_{\pm}$ is
unidentifiable from EIS data and thus only estimates from NLEIS data are
shown. Also, although NLEIS data can be used to estimate the typical non-
dimensional charge transfer resistance $\chi_{\pm}$, EIS data cannot, so to
compare methods we plot the charge transfer resistances $R_{\pm}^{(0)}$
instead.

Figs. 8 and 9 show that the mean and variance of the kinetic parameter
estimators $R_{\pm}^{(0)},C_{\pm}$ are similar in all cases. The NLEIS-based
estimator of $\beta_{+}$ has low variance, but that for $\beta_{-}$ has high
variance. This is due to the resistance of the negative electrode
($R_{-}^{(0)}$) being much lower than that of the positive, hence its kinetic
cardioid spiral in $Z_{2}^{*}$ is much smaller in magnitude, and affected more
by noise—see Fig. 8$(j)$-$(k)$. The discrepancy between the means and the true
values (i.e. estimator bias) is due to the small-current-amplitude
approximation we used to derive the model impedances. We determined this by
fitting to synthetic data generated using smaller current amplitudes and found
that the bias reduced. This bias thus represents model error rather than
noise, as the data generation used a finite amplitude.

The diffusion timescales $\tau_{\mathrm{d},\pm}$ are more influenced by noise,
with the NLEIS estimates showing larger variance than the EIS ones. However,
this variance is comparable to the model error—the EIS estimates of
$\tau_{\mathrm{d},\pm}$ can be overly confident (low variance), despite having
a similar model error to NLEIS (Fig. 9 for 70% DoD).

Figure 8: (50% DoD) Parameter distributions (histograms) for
$\hat{\bm{\theta}}_{\mathrm{NLEIS}}$ from repeated fitting to synthetic NLEIS
data, showing sensitivity to error $(a)-(h)$. Distributions for fitting
$\hat{\bm{\theta}}_{\mathrm{EIS}}$ to just EIS data (fundamental only) also
shown. Dashed lines are mean values, with true ones shown in green. Estimates
for $\beta_{\pm}$ not shown for EIS as they are not identifiable. Example
impedance fit shown in $(i)$-$(k)$, including blown-up plot $(k)$ of $Z_{2}$
near the origin.
Figure 9: (70% DoD) Parameter distributions (histograms) for
$\hat{\bm{\theta}}_{\mathrm{NLEIS}}$ from repeated fitting to synthetic NLEIS
data, showing sensitivity to error. See caption for Fig. 8.

#### Log-likelihood analysis of diffusion timescales

A clearer understanding of the practical identifiability of diffusion time
$\tau_{\mathrm{d},\pm}$ can be reached by looking at the log-likelihoods (i.e.
the cost functions associated with the parameter estimation problem) for each
harmonic directly, $l_{1}(\bm{\theta})$, $l_{2}(\bm{\theta})$ in the
$(\tau_{\mathrm{d},+},\tau_{\mathrm{d},-})$ plane. Fig. 10 shows these (and
the full log-likelihood $l_{12}=l_{1}+l_{2}$) at 3 DoDs, with kinetic
parameters fixed at their true values.

The shape of the $l_{1}$ landscape (Fig. 10$(a,d,g)$) depends on OCP slopes
$U_{\pm}^{(0)\prime}$, as discussed in detail by Bizeray et al. [21] (for
Fickian diffusion, but the structure is qualitatively the same). When one of
the OCP slopes $U_{\pm}^{(0)\prime}$ is small, a valley in $l_{1}$ appears
through the minimum, showing poor identifiability of the corresponding
$\tau_{\mathrm{d},\pm}$. Of the results shown, the negative electrode OCP
slope $|U_{-}^{(0)\prime}|$ is largest ($\approx 17.9$) at 50% DoD, where
$\tau_{\mathrm{d},-}$ identifiability is best, and smallest ($\approx 0.36$)
at 30% DoD, where it is worst. But the positive electrode OCP slope
$|U_{+}^{(0)\prime}|$ is largest ($\approx 18.0$) at 30% DoD and decreases at
50% DoD ($\approx 6.5$) and 70% DoD ($\approx 2.2$), decreasing
$\tau_{\mathrm{d},+}$ identifiability. At 70% DoD,
$|U_{-}^{(0)\prime}|\approx|U_{+}^{(0)\prime}|$, giving electrode symmetry
(see the Structural Parameter Identifiability section), and 2 local minima in
$l_{1}$. The red line on the $l_{1}$ plots is the relationship between
$\tau_{\mathrm{d},-}$ and $\tau_{\mathrm{d},+}$ if both slopes
$U_{-}^{(0)\prime},U_{+}^{(0)\prime}\to 0$, in which case (86) reduces to, at
leading order,

$\displaystyle
A_{+}\sqrt{\tau_{\mathrm{d},+}}+A_{-}\sqrt{\tau_{\mathrm{d},-}}$
$\displaystyle=\mbox{const.}$ $\displaystyle\mbox{where }A_{\pm}$
$\displaystyle=\xi_{\pm}\sqrt{\frac{|U_{\pm}^{(0)\prime}|}{c_{\pm}^{(0)}}}.$
(89)

This gives a visual guide to the predominant structure of the $l_{1}$
landscape.

For the log-likelihood $l_{2}(\bm{\theta})$ of the second harmonic (Fig.
10$(b,e,h)$), its landscape depends on
$U_{\pm}^{(0)\prime},U_{\pm}^{(0)\prime\prime}$ and $R_{\pm}^{(0)\prime}$ and
so is more difficult to predict from the formulae, but crucially it is
distinct from the landscape of $l_{1}$. Thus, when it is combined with $l_{1}$
to produce the total log-likehood $l_{12}$, the identifiability of the pair
$(\tau_{\mathrm{d},+},\tau_{\mathrm{d},-})$ is improved. For example, at 50%
and 70% DoD, identifiability of $\tau_{\mathrm{d},+}$ is improved, and the
second minimum (at 70%) eliminated. In each case, including $l_{2}$ moves the
global minimum closer to the true value. We make an additional remark that the
minima of $l_{1}$ in Fig. 10 are further from the true value than the
estimates in Fig. 8. This is because, in the latter, the kinetic parameters
were fitted simultaneously which allowed a better estimate of the diffusion
times (i.e. reduced the bias due to the small current amplitude expansion). As
kinetic parameters were fixed at their true values in Fig. 10, this bias is
increased for the diffusion times. Using fitted kinetic parameters in Fig. 10
instead, the minima move closer to the true values. However, this highlights
the high sensitivity of linear EIS estimates to errors in the kinetic
parameters.

Lastly, relevant for fitting to experimental data in the next section, we
consider the case where OCP curvature $U_{\pm}^{(0)\prime\prime}$ is not
known, due to measurement error and difficulties associated with approximating
second derivatives numerically from noisy data. In this case, Fig. 11 shows
$l_{1},l_{2},l_{12}$ where, for each pair of values
$(\tau_{\mathrm{d},+},\tau_{\mathrm{d},-})$, we have minimised over
$U_{-}^{(0)\prime\prime}$ and $U_{+}^{(0)\prime\prime}$—this corresponds to
constructing “profile likelihoods” where
$U_{-}^{(0)\prime\prime},U_{+}^{(0)\prime\prime}$ are nuisance parameters that
are not of interest [30, 31]. Note that $l_{1}$ is unchanged from Fig. 10 as
it does not depend on $U_{\pm}^{(0)\prime\prime}$. Even here, with no
knowledge of $U_{\pm}^{(0)\prime\prime}$, the second harmonic provides enough
information to improve identifiability over the fundamental alone.

In Fig. 11, the dashed overlaid line corresponds to
$\tau_{\mathrm{d},-}=(c_{+}^{(0)}U_{+}^{(0)\prime})\tau_{\mathrm{d},+}/(c_{-}^{(0)}U_{-}^{(0)\prime})$
(or $D_{-}^{(0)}=D_{+}^{(0)}$) along which the pair
$U_{-}^{(0)\prime\prime},U_{+}^{(0)\prime\prime}$ are unidentifiable. This was
verified by investigating the
$(U_{-}^{(0)\prime\prime},U_{+}^{(0)\prime\prime})$ plane when
$(\tau_{\mathrm{d},+},\tau_{\mathrm{d},-})$ are on this dashed line. However,
this only introduces small numerical irregularities which are irrelevant if
far from the global minimum.

Figure 10: Synthetic data: Log-likelihood plotted in the $\tau_{d,+}$ and
$\tau_{d,-}$ plane (other parameters fixed at true values). DoD (indicated)
increases from top row to bottom row. Left column is log-likelihood $l_{1}$
for $Z_{1}^{(1)}$ only; middle column is log-likelihood $l_{2}$ for
$Z_{2}^{(2)}$ only; right column is the total log-likelihood
$l_{12}=l_{1}+l_{2}$. For 30%, 50%, 70% DoD we have $-U_{-}^{(0)\prime}\approx
0.4,17.9,1.6$, and $-U_{+}^{(0)\prime}\approx 18.0,6.5,2.2$. Red line is (89).
Figure 11: Synthetic data: Log-likelihood plotted in the $\tau_{d,+}$ and
$\tau_{d,-}$ plane. Same as Fig. 10, but at each ($\tau_{d,+},\tau_{d,-}$)
pair, we minimize $l_{2}$ over $U_{+}^{(0)\prime\prime}$ and
$U_{-}^{(0)\prime\prime}$. Dashed line is
$\tau_{\mathrm{d},-}=(c_{+}^{(0)}U_{+}^{(0)\prime})\tau_{\mathrm{d},+}/(c_{-}^{(0)}U_{-}^{(0)\prime})$
(or $D_{-}=D_{+}$) where $U_{\pm}^{(0)\prime\prime}$ are unidentifiable.

###  5.4 Parameter estimation from experimental data

We have explored the parameter estimation approach and advantages of NLEIS
using synthetic data, and now turn to experimental data collected from an NMC
cathode/graphite anode cell—see the Experimental Methods section for details.
We consider data collected at 10%, 20%, …, 90% DoD, at each of which we
consider EIS ($Z_{1}^{*}$ only) and NLEIS ($Z_{1}^{*}$ and $Z_{2}^{*}$) model
fits.

The fits to EIS data are shown first, in Fig. 12, with the associated
parameter estimates in Fig. 15. The model is able to fit the semi-circular
reaction kinetics well, although we do not expect to capture the very high
frequency behaviour as our model does not include phenomena such as transport
through SEI layers, relevant for $\omega^{*}\gtrsim 1$ kHz. We can extract the
series resistance as the high frequency intercept with the real axis, giving
$R_{s}^{*}=$ $9.2\text{\,}\mathrm{m\SIUnitSymbolOhm}$ (and dimensionless
$R_{s}=0.358$). As discussed in the Structural Parameter Identifiability
section, estimating kinetics from EIS data suffers from an inability to
determine which electrode the parameters correspond to. Hence, we made the
assumption that the larger charge transfer resistance (radius of semi-circle)
was due to the positive electrode, but this was arbitrary and did not affect
the goodness of fit. Note also that $\beta_{-}$ and $\beta_{+}$ cannot be
determined at all, so we set them equal to 1/2, as is typically done
implicitly in other studies.

The diffusion tails are fitted reasonably well. However the corresponding
timescales show extreme variation with DoD, up to 4 orders of magnitude—as is
usually seen when estimating diffusion parameters from EIS or GITT data [32,
33, 34]. As we will see, NLEIS can be a more sensitive tool for probing the
validity of the underlying diffusion model.

Next, fits to NLEIS data ($Z_{1}^{*}$ and $Z_{2}^{*}$ simultaneously) are
shown in Fig. 13 (with parameters in Fig. 15). We consider two cases: (i)
using values of $U_{\pm}^{(0)\prime\prime}$ calculated from Fig. 5, and; (ii)
fitting $U_{\pm}^{(0)\prime\prime}$, assuming they are not known _a priori_ ,
or are inaccurate. In both cases, the reaction kinetics in each harmonic are
captured very well for each DoD, as shown by the semi-circles in $Z_{1}^{*}$
and cardioid spirals in $Z_{2}^{*}$. The resulting values for
$\beta_{\pm},C_{\pm},\chi_{\pm}$ also do not vary appreciably across DoD
except $\chi_{+}$, suggesting the functional form of the exchange current
density (6) is less appropriate for that electrode. The charge transfer
coefficients $\beta_{\pm}$ both only slightly deviate from the symmetric value
of 1/2 for these cells. Here $\beta_{+}$ deviates most from 1/2 and is
responsible for the predominant cardioid in $Z_{2}^{*}$. The cardioid from the
negative electrode is small and only visible by zooming in closer to the
origin—see Fig. 14.

Regarding the low frequency diffusion tails, the model fits are less
convincing. When $U_{\pm}^{(0)\prime\prime}$ is assumed known, the tail in
$Z_{2}^{*}$ may point in the opposite direction (10%, 20%, 50% DoD) or have
the wrong length (30%, 40%). Also, the tail in $Z_{1}^{*}$ is now too short in
some cases (40%, 50%, 60%). However, if OCP curvatures
$U_{\pm}^{(0)\prime\prime}$ are allowed to be fitted instead, the fit to the
tail in $Z_{2}^{*}$ is significantly improved. Also, the tail in $Z_{1}^{*}$
and the fit to the kinetic portions are improved; see the zoom-ins of Fig. 14.
This suggests, even if only the _kinetic_ parameters are of interest, fitting
$U_{\pm}^{(0)\prime\prime}$ may be better than prescribing inaccurate values.
The resulting diffusion timescales (Fig. 15$(c,d)$ also show the least
variation with DoD for NLEIS data fits when $U_{\pm}^{(0)\prime\prime}$ is
fitted (red), an improvement over the EIS data fits (blue), and hence closer
to the model assumption that it is does not vary with DoD.

Figure 12: Fits of the model fundamental impedance $Z_{1}^{(1)*}$ to
corresponding experimental fundamental impedance data. Units are in Ohms
(VA-1). Fits performed at each DoD in the range 10%, 20%,… 90% are shown, for
the frequency range $f^{*}=0.01-500$ Hz. Parameter estimates are given in Fig.
15. Figure 13: Simultaneous fits of the model impedances $Z_{1}^{(1)*}$
(fundamental) and $Z_{2}^{(2)*}$ (second harmonic) to corresponding
experimental data. $Z_{1}^{*},Z_{2}^{*}$ units are in VA-1 (Ohms) and VA-2.
Two cases, where $U_{\pm}^{\prime\prime}$ is either fixed (assumed known) or
chosen as a fitting parameter are shown. Fits performed at each DoD in the
range 10%, 20%,… 90% separately, increasing from top left to the bottom right
panel. Parameter estimates are given in Fig. 15. Figure 14: Successive zoom-
ins on selected $Z_{2}^{*}$ impedances in Fig. 13. The top row is at 10% DoD,
and the bottom row is at 60% DoD. Figure 15: Parameter values from fitting to
experimental EIS data ($Z_{1}^{*}$ only) or NLEIS data ($Z_{1}^{*}$ and
$Z_{2}^{*}$) at different DoDs, from 10% to 90%. Kinetic parameters are in the
left two columns ($a,b,e,f,i,j$), with diffusion timescales and OCP parameters
in the remaining two; $\beta_{+}$ and $\beta_{-}$ were taken as 0.5 in the EIS
fits.

#### Log-likelihood analysis of diffusion timescales

The fits of diffusion times to experimental data can be further probed by
examining the underlying likelihood landscape in the
$(\tau_{\mathrm{d},+},\tau_{\mathrm{d},-})$ plane. As we did for synthetic
data, we plot here the log-likelihoods in Fig. 16, minimising over
$U_{\pm}^{(0)\prime\prime}$. The DoDs chosen (50% and 60%) are ones at which
the timescale estimates from EIS vs. NLEIS differ significantly. The rest of
the parameters are fixed at their estimated values from Fig. 15.

The $l_{1}$ landscape ($Z_{1}^{*}$) clearly shows unidentifiability given by
the red line, i.e. (89). Here this is likely due to insufficient data at lower
frequencies, as (89) also holds if only the high frequency approximation (48)
is taken in the Warburg impedance. In the experiments, we were limited to a
minimum allowable frequency of $0.01\text{\,}\mathrm{H}\mathrm{z}$ by the
Kikusui PBZ60-6.7 power supply. This was not an issue when estimating from
synthetic data (Fig. 11) as we were able to generate data at arbitrarily low
frequencies—observe that the diffusion tail in Fig. 8($i$) becomes close to
vertical at the lowest frequencies. This unidentifiability valley was
therefore not present, with one (or 2) local minima present instead.

The $l_{2}$ landscape ($Z_{2}^{*}$) in Fig. 16 has its minima considerably far
from where $l_{1}$ predicts they should be. Further, the depth of $l_{2}$
minima overwhelms those in $l_{1}$ when combined in $l_{12}=l_{1}+l_{2}$. This
differs significantly from the behaviour for synthetic data where the model
exactly matches that used to generate the data, i.e. there is no model error,
only approximation error due to the expansion in $\hat{I}$, and noise. In that
case, the results showed that the fundamental is more dominant, with deeper
minima in $l_{1}$ than $l_{2}$.

This inconsistency between the minima of $l_{1}$ and $l_{2}$ suggests that
there may be deficiencies in the diffusion model even though the model is able
to to fit the data very well. Here we considered a nonlinear diffusion model,
but to alleviate these issues one may need to consider other approaches, e.g.
those incorporating phase-field dynamics (such as Cahn–Hilliard) which have
shown to better represent lithiation in some electrode materials, e.g.
graphite [35]. On the other hand, it may be that a single representative
particle for each electrode is inadequate, or that large-scale geometric
effects are important, or the inclusion of transport in the electrolyte is
necessary. These results suggest NLEIS could be a promising tool for model
selection. Other dynamics that may also relate to this inconsistency are those
in the electrolyte, e.g. diffusion and migration of lithium, which have been
neglected here.

Figure 16: Log-likelihoods in the diffusion timescale plane
$\tau_{\mathrm{d},+}$-$\tau_{\mathrm{d},-}$, using experimental data. Same as
Fig. 10 but now at each point we have minimised over $U_{-}^{\prime\prime}$
and $U_{+}^{\prime\prime}$.

###  5.5 Model validation in the time domain

We have parameterized the nonlinear SPM using experimental NLEIS (frequency-
domain) data and now proceed to validate the model performance and parameters
by numerical solution in the time domain. We considered a current demand based
on an electric vehicle drive cycle (see [36]) with measurements taken on the
same model of Kokam cell. Current was specified at 1 second intervals, with
peaks reaching $4.4\text{\,}\mathrm{A}$ (or C-rates of 6), and initiated at 0%
DoD. The model (15)-(24) was solved using this profile as an input, which was
interpolated using piecewise-cubic Hermite interpolating polynomials (PCHIP)
to ensure no overshooting of the data. Further details of numerical methods
are given in the Methods.

We compared the model performance using two parameter sets: estimates from EIS
data, and estimates from NLEIS data (see Fig. 15). We used NLEIS parameters
for the case where $U_{\pm}^{\prime\prime}$ were considered unknown. For
$\beta_{\pm},C_{\pm},\tau_{\mathrm{d},\pm}$, we chose single representative
values by averaging over the full range of DoDs (with the average of
$\tau_{\mathrm{d},\pm}$ taken on $\log\tau_{\mathrm{d},\pm}$). For
$\chi_{\pm}$, we converted Fig. 15 data from DoD to stoichiometry $c_{\pm}$
and interpolated (PCHIP) to give the functions
$\chi_{+}(c_{+}),\chi_{-}(c_{-})$.222$\chi_{\pm}$ were originally assumed
constant in the model, but making them depend on $c_{\pm}$ accounts for errors
in the functional form of $R_{\pm}(c_{\pm})$, (24). The parameter values used
are summarized in Table 1.

Comparison between the model and experiment is shown in Fig. 17. In general,
the model results with either parameter set (EIS or NLEIS) agree excellently
with the measured data over the full duration ($\approx 20$ min), with no
parameter tuning performed.333The series resistance $R_{s}^{*}$ was changed
before simulation to 20 m$\Omega$ (taken from Fig. 4 of [21])) since the drive
cycle measurements were performed using a different test rig, but no
subsequent changes to $R_{s}^{*}$ were made. Even the EIS predictions show an
RMSE of only $15.2\text{\,}\mathrm{m}\mathrm{V}$ and no overall drift in error
as time increases. This is a marked improvement over the linear SPM of Bizeray
_et al._ [21] where, for the same cells and validation data, they required the
linearized OCPs to be adjusted to correct for drift in the DoD . Our model
exhibited no discernable drift and we considered a longer drive cycle
simulation than Bizeray _et al._ [21] ($\approx 17$ mins versus $9$ mins).

The predictions from NLEIS, however, agree even better with the measured
voltage data, with a lower RMSE of $11.1\text{\,}\mathrm{m}\mathrm{V}$. We
conclude that NLEIS can provide better estimates of the reaction kinetic
parameters ($\chi_{\pm}$, $\beta_{\pm}$, and $C_{\pm}$). Looking at the
stoichiometries on the surface of electrode particles (Fig. 17$(d)$) we see
appreciable differences only in the negative electrode, but this has minimal
effect on the voltage since $U_{-}(c_{-})$ is mostly flat over the range of
$c_{-}$ operated in. Hence, diffusion timescale estimates from either method
are sufficient for good voltage accuracy in the time domain, but predictions
of internal states may differ greatly, which will have consequences for
degradation modelling.

Figure 17: Model validation in the time domain, corresponding to a drive cycle of prescribed current $(a)$. The terminal voltage from experiment [36] and two model simulations, using parameter values (Table 1) from fitting to EIS or NLEIS data, are in $(b)$. Also shown are $(c)$ the absolute voltage error, and $(d)$ surface stoichiometries $c_{s,\pm}$ in each electrode for each simulation. Parameter groups | Negative (-) | Positive (+)
---|---|---
EIS estimate | NLEIS estimate | EIS estimate | NLEIS estimate
$\chi_{\pm}$ |  Fig. 15$(a)$ |  Fig. 15$(a)$ |  Fig. 15$(b)$ |  Fig. 15$(b)$
$\beta_{\pm}$ | 0.5 | 0.496 | 0.5 | 0.482
$C_{\pm}$ | $1.528\times 10^{-3}$ | $2.589\times 10^{-3}$ | $1.747\times 10^{-2}$ | $2.645\times 10^{-2}$
$\tau_{\mathrm{d},\pm}$ | $1.571\times 10^{4}$ | $1.570\times 10^{3}$ | $1.227\times 10^{5}$ | $1.185\times 10^{5}$
$\xi_{\pm}$ | $8.126\times 10^{-5}$ | $7.414\times 10^{-5}$
$c_{\pm}^{0\%}$ | 0.664 | 0.40
Table 1: Parameter group values estimated from EIS and NLEIS data, which are
used for model validation in Fig. 17. The $\chi_{\pm}$ parameters are
interpolated on the values in Fig. 15$(a,b)$ for different DoDs (NLEIS
estimates used are the red markers). From our impedance measurements the
series resistance is $R_{s}^{*}=9.2$ mV ($R_{s}=0.358$), but for the model
validation (Fig. 17) we use $R_{s}^{*}=20$ mV and $R_{s}=0.779$.

##  6 Conclusions

In this paper, we presented a comprehensive study on the application of NLEIS
to a standard physical model of a lithium-ion battery: the single particle
model. We derived a suite of analytical formulae for the fundamental harmonic
and second harmonic impedances resulting from this model, and then interpreted
each contribution to the formulae physically. Signatures of various processes
and model nonlinearities were identified in the higher harmonics, in
particular: (i) asymmetric charge transfer in Butler–Volmer kinetics; (ii)
concentration dependence of the exchange current density; (iii) nonlinearity
of the electrode OCPs; and (iv) nonlinearity of solid-state diffusion (i.e.
concentration dependence of the diffusivity). As a result, we demonstrated how
the Nyquist plot of the second harmonic impedance may be interpreted in a
similar manner to that of the fundamental, but built from different features
and therefore containing additional information about the system and its
parameters. This information was then exploited for parameter estimation from
experimental data that was collected from commercial cells.

From our nonlinear impedance formulae, we explored how the second harmonic
impedance in combination with the usual fundamental impedance can be used to
improve parameter identifiability, both structurally and practically, over the
use of the fundamental alone (EIS). Structurally, NLEIS can identify new
parameters not possible from EIS, e.g. charge transfer coefficients (usually
assumed to be $1/2$ for Li-ion batteries), and improve the identifiability of
the remaining parameters. In particular, positive-negative electrode
symmetries are eliminated, allowing all parameters to be associated with the
correct electrode, and the curvature (if nonzero) of the OCP can help identify
a diffusion timescale which would otherwise be unidentifiable from EIS (e.g.,
due to a “flat” OCP change with respect to stoichiometry).

To investigate the impact of NLEIS on practical identifiability, we
implemented maximum likelihood parameter estimation. First, we demonstrated
the procedure using synthetic data and bootstrap resampling, showing that one
can robustly determine all the model parameter groups (kinetic and diffusive)
from NLEIS data. The second harmonic has a lower signal-to-noise ratio
resulting in larger variances for some parameters, but these uncertainties are
comparable in magnitude to model error, e.g. due to the small excitation
amplitude approximation.

Finally, our parameter estimation procedure was applied to experimental NLEIS
data that was collected from a commercial cell. We could fit the impedance
data very well, capturing many of the key features of both harmonics. All the
kinetic parameters were determined, and the charge transfer coefficients were
found to deviate slightly from $1/2$, albeit remaining within 10%. However,
the diffusion timescale estimates from EIS and NLEIS data differed greatly and
the latter, via analysis of the likelihood landscape, suggested a discrepancy
in the diffusion model—a discrepancy otherwise not revealed from EIS data
alone. Nonetheless, validation of the model and its parameterization from
either EIS or NLEIS data was conducted in the time domain with a drive cycle
experiment, where an improvement in accuracy was observed for the NLEIS
parameterization, showing that the technique is beneficial when the model is
used in practice.

Areas of further work could consider the impacts of several physical
mechanisms that were neglected here, such as transport through SEI layers and
in the electrolyte. To keep a similar level of model complexity, one could
consider the latter by extending the current analysis to the single particle
model with electrolyte (SPMe) [26]. The electrolyte introduces additional
diffusive effects which may address the difficulties observed here with
fitting the diffusion tails in a consistent manner. Other solid-state
transport models, such as phase-field models [35], could also be considered,
and NLEIS used as a model selection tool. In addition, extending the minimum
frequency measured experimentally could improve diffusion time estimates.
Finally, metrics to quantify model parsimony could be considered, i.e. to
analyse whether there are too many (or too few) parameters, and to compare
different model structures.

##  7 Acknowledgements

This research was funded in whole, or in part, by the EPSRC Faraday
Institution Multiscale Modelling project (EP/S003053/1, grant number FIRG025).
The authors also wish to thank Robert Timms and Steven Psaltis for helpful
discussions. For the purpose of Open Access, the authors have applied a CC BY
public copyright licence to any Author Accepted Manuscript (AAM) version
arising from this submission.

## References

  * [1] C. Iclodean, B. Varga, N. Burnete, D. Cimerdean, B. Jurchiş, Comparison of different battery types for electric vehicles, IOP Conference Series: Materials Science and Engineering 252 (2017) 012058.  doi:10.1088/1757-899x/252/1/012058.
  * [2] M. Doyle, T. M. Fuller, J. Newman, Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell., Journal of the Electrochemical Society 140 (6) (1993) 1526–1533.
  * [3] Y. Barsukov, J. R. Macdonald, Electrochemical Impedance Spectroscopy, American Cancer Society, 2012, pp. 1–17.  doi:10.1002/0471266965.com124.
  * [4] F. Ciucci, Modeling electrochemical impedance spectroscopy, Current Opinion in Electrochemistry 13 (2019) 132 – 139, fundamental and Theoretical Electrochemistry: Physical and Nanoelectrochemistry.  doi:https://doi.org/10.1016/j.coelec.2018.12.003.
  * [5] M. D. Murbach, D. T. Schwartz, Extending newman pseudo-two-dimensional lithium-ion battery impedance simulation approach to include the nonlinear harmonic response, Journal of The Electrochemical Society 164 (11) (2017) E3311–E3320.  doi:10.1149/2.0301711jes.
  * [6] N. Harting, N. Wolff, F. Röder, U. Krewer, Nonlinear frequency response analysis (nfra) of lithium-ion batteries, Electrochimica Acta 248 (2017) 133–139.  doi:https://doi.org/10.1016/j.electacta.2017.04.037.
  * [7] T. J. McDonald, S. Adler, (invited) theory and application of nonlinear electrochemical impedance spectroscopy, ECS Transactions 45 (1) (2012) 429–439.  arXiv:http://ecst.ecsdl.org/content/45/1/429.full.pdf+html, doi:10.1149/1.3701334.
URL http://ecst.ecsdl.org/content/45/1/429.abstract

  * [8] J. Wilson, D. Schwartz, S. Adler, Nonlinear electrochemical impedance spectroscopy for solid oxide fuel cell cathode materials, Electrochimica Acta 51 (8) (2006) 1389 – 1402, electrochemical Impedance Spectroscopy.  doi:https://doi.org/10.1016/j.electacta.2005.02.109.
  * [9] J. R. Wilson, M. Sase, T. Kawada, S. B. Adler, Measurement of Oxygen Exchange Kinetics on Thin-Film La0.6Sr0.4CoO3 Using Nonlinear Electrochemical Impedance Spectroscopy, Electrochemical and Solid-State Letters 10 (5) (2007) B81–B86.  doi:10.1149/1.2710178.
  * [10] T. Kadyk, R. Hanke-Rauschenbach, K. Sundmacher, Nonlinear frequency response analysis of pem fuel cells for diagnosis of dehydration, flooding and co-poisoning, Journal of Electroanalytical Chemistry 630 (1) (2009) 19 – 27.  doi:https://doi.org/10.1016/j.jelechem.2009.02.001.
  * [11] T. Kadyk, R. Hanke-Rauschenbach, K. Sundmacher, Nonlinear frequency response analysis for the diagnosis of carbon monoxide poisoning in pem fuel cell anodes, Journal of Applied Electrochemistry 41 (9) (2011) 1021–1032.  doi:10.1007/s10800-011-0298-8.
  * [12] N. Xu, J. Riley, Nonlinear analysis of a classical system: The double-layer capacitor, Electrochemistry Communications 13 (10) (2011) 1077 – 1081.  doi:https://doi.org/10.1016/j.elecom.2011.07.003.
  * [13] N. Xu, D. J. Riley, Nonlinear analysis of a classical system: The faradaic process, Electrochimica Acta 94 (2013) 206 – 213.  doi:https://doi.org/10.1016/j.electacta.2013.01.141.
  * [14] F. Fasmin, R. Srinivasan, Review: Nonlinear electrochemical impedance spectroscopy, Journal of The Electrochemical Society 164 (7) (2017) H443–H455.  doi:10.1149/2.0391707jes.
  * [15] 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–A2765.  doi:10.1149/2.0711811jes.
  * [16] N. Harting, N. Wolff, U. Krewer, Identification of lithium plating in lithium-ion batteries using nonlinear frequency response analysis (nfra), Electrochimica Acta 281 (2018) 378–385.  doi:https://doi.org/10.1016/j.electacta.2018.05.139.
  * [17] N. Wolff, N. Harting, M. Heinrich, F. Röder, U. Krewer, Nonlinear frequency response analysis on lithium-ion batteries: A model-based assessment, Electrochimica Acta 260 (2018) 614–622.  doi:https://doi.org/10.1016/j.electacta.2017.12.097.
  * [18] N. Wolff, N. Harting, M. Heinrich, U. Krewer, Nonlinear frequency response analysis on lithium-ion batteries: Process identification and differences between transient and steady-state behavior, Electrochimica Acta 298 (2019) 788–798.  doi:https://doi.org/10.1016/j.electacta.2018.12.107.
  * [19] C. Fan, K. O’Regan, L. Li, E. Kendrick, W. Widanage, Frequency domain non-linear characterization and analysis of lithium-ion battery electrodes, Journal of Energy Storage 36 (2021) 102371.  doi:https://doi.org/10.1016/j.est.2021.102371.
  * [20] C. Fan, T. R. Grandjean, K. O’Regan, E. Kendrick, W. D. Widanage, Understanding non-linearity in electrochemical systems using multisine-based non-linear characterization, Transactions of the Institute of Measurement and Control (2021).  doi:10.1177/01423312211045991.
  * [21] A. M. Bizeray, J. Kim, S. R. Duncan, D. A. Howey, Identifiability and parameter estimation of the single particle lithium-ion battery model, IEEE Transactions on Control Systems Technology 27 (5) (2019) 1862–1877.  doi:10.1109/TCST.2018.2838097.
  * [22] S. Atlung, K. West, T. Jacobsen, Dynamic aspects of solid solution cathodes for electrochemical power sources, Journal of The Electrochemical Society 126 (8) (1979) 1311–1321.  doi:10.1149/1.2129269.
  * [23] G. Ning, B. N. Popov, Cycle life modeling of lithium-ion batteries, Journal of The Electrochemical Society 151 (10) (2004) A1584–A1591.  doi:10.1149/1.1787631.
  * [24] J. S. Horner, G. Whang, D. S. Ashby, I. V. Kolesnichenko, T. N. Lambert, B. S. Dunn, A. A. Talin, S. A. Roberts, Electrochemical modeling of gitt measurements for improved solid-state diffusion coefficient evaluation, ACS Applied Energy Materials 4 (10) (2021) 11460–11469.  doi:10.1021/acsaem.1c02218.
  * [25] E. McTurk, C. R. Birkl, M. R. Roberts, D. A. Howey, P. G. Bruce, Minimally invasive insertion of reference electrodes into commercial lithium-ion pouch cells, ECS Electrochemistry Letters 4 (12) (2015) A145–A147.  doi:10.1149/2.0081512eel.
  * [26] S. G. Marquis, V. Sulzer, R. Timms, C. P. Please, S. J. Chapman, An asymptotic derivation of a single particle model with electrolyte, Journal of The Electrochemical Society 166 (15) (2019) A3693.  doi:10.1149/2.0341915jes.
  * [27] Kokam Cell SLPB533459H4 Data Sheet (Accessed 2022).
URL https://www.custompower.com/documents/SLPB533459H4.pdf

  * [28] R. Bellman, K. Åström, On structural identifiability, Mathematical Biosciences 7 (3) (1970) 329–339.  doi:https://doi.org/10.1016/0025-5564(70)90132-X.
  * [29] B. Efron, R. Tibshirani, An Introduction to the Bootstrap, Chapman and Hall/CRC, 1993.
  * [30] D. Cox, Principles of Statistical Inference, Cambridge University Press, 2006.
  * [31] M. J. Simpson, R. E. Baker, S. T. Vittadello, O. J. Maclaren, Practical parameter identifiability for spatio-temporal models of cell invasion, Journal of The Royal Society Interface 17 (164) (2020) 20200055.  doi:10.1098/rsif.2020.0055.
  * [32] M. Ecker, S. Käbitz, I. Laresgoiti, D. U. Sauer, Parameterization of a physico-chemical model of a lithium-ion battery: II. Model validation, Journal of The Electrochemical Society 162 (9) (2015) A1849–A1857.  doi:10.1149/2.0541509jes.
  * [33] J. Schmalstieg, C. Rahe, M. Ecker, D. U. Sauer, Full cell parameterization of a high-power lithium-ion battery for a physico-chemical model: Part i. physical and electrochemical parameters, Journal of The Electrochemical Society 165 (16) (2018) A3799–A3810.  doi:10.1149/2.0321816jes.
  * [34] 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.  doi:10.1149/1945-7111/ab9050.
  * [35] T. R. Ferguson, M. Z. Bazant, Phase transformation dynamics in porous battery electrodes, Electrochimica Acta 146 (2014) 89 – 97.  doi:https://doi.org/10.1016/j.electacta.2014.08.083.
  * [36] S. Zhao, S. R. Duncan, D. A. Howey, Observability analysis and state estimation of lithium-ion batteries in the presence of sensor biases, IEEE Transactions on Control Systems Technology 25 (1) (2017) 326–333.  doi:10.1109/TCST.2016.2542115.