Second-Harmonic Nonlinear Electrochemical Impedance Spectroscopy: Part I. Analytical Theory and Equivalent Circuit Representations for Planar and Porous Electrodes

Electrochemical impedance spectroscopy (EIS) experiments impose a single frequency alternating current excitation to the electrode under study, represented here as

$$\begin{eqnarray}&&I=\unicode{x02206}I\exp \left(j\omega t\right)\end{eqnarray} \tag{ 5 }$$

where $j=\sqrt{-1},$ $\unicode{x02206}I$ is the amplitude of the current modulation (Amp) with frequency $\omega$ (radians/s) For sufficiently small excitation amplitude $\unicode{x02206}I$ , the voltage output for a stationary system is dominated by the linear response of the interface and produces a modulated voltage at the excitation frequency $\omega .$ Murbach et al. has shown theoretically, ²³ and validated experimentally, ²² that moderately larger amplitude single-frequency current excitations can generate a weakly nonlinear modulated voltage response of the form

$$\begin{eqnarray}&&E={\widetilde{E}}_{1}\exp \left(j\omega t\right)+{\widetilde{E}}_{2}\exp \left(j2\omega t\right)+h.o.t.\end{eqnarray} \tag{ 6 }$$

where the voltage output is at the excitation frequency $\omega$ and the second harmonic $2\omega .$ To be consistent with Eq. 1, we ignore higher order terms generated at increasingly large current modulation amplitudes, but can direct the interested reader elsewhere for the formal solution structure for h.o.t. ^12,26,27 As we have shown previously, ²² the linear and second harmonic impedances from the voltage spectrum represented by Eq. 6 are

$$\begin{eqnarray}&&{\widetilde{Z}}_{1}\left(\omega \right)=\frac{{\widetilde{E}}_{1}}{\unicode{x02206}I}\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{\widetilde{Z}}_{2}\left(\omega \right)=\frac{{\widetilde{E}}_{2}}{{(\unicode{x02206}I)}^{2}}\end{eqnarray} \tag{ 8 }$$

where ${\widetilde{Z}}_{1}$ is the complex frequency-dependent (linear) EIS and ${\widetilde{Z}}_{2}$ is the complex frequency-dependent second harmonic nonlinear EIS (2nd-NLEIS) function.

Mathematically, the EIS function ${\widetilde{Z}}_{1}\left(\omega \right)$ computed in Eq. 7 is identical to the comparable (linear) EIS spectrum when Eq. 6 represents the voltage spectrum; practically, we have shown that achieving Eq. 6 experimentally is straightforward. ²² In contrast, when the current amplitude is large enough to generate half-cell voltage spectra with three, four, or higher multiples of the excitation frequency, then ${\widetilde{Z}}_{1}\left(\omega \right)$ computed in Eq. 7 will have additional nonlinear contributions that deviate from the traditional (linear) response. Here, we only consider moderate modulation amplitudes consistent with Eq. 6.

A planar electrochemical interface with time-invariant and uniform concentrations throughout the solid and electrolyte phases is governed by Eq. 3 with constant ${E}_{n}.$ For simplicity, we set ${E}_{n}=0$ and re-write Eq. 3 in terms of physically interpretable parameters common to a Randles circuit analysis, namely,

$$\begin{eqnarray}&&I=\frac{E}{{R}_{{ct}}}+\frac{\varepsilon f{E}^{2}}{{R}_{{ct}}}+{C}_{{dl}}\frac{{dE}}{{dt}}\end{eqnarray} \tag{ 9 }$$

where ${R}_{{ct}}=\frac{1}{{{Ai}}_{0}f}$ is the interfacial charge transfer resistance (${\rm{\Omega }}$), and ${C}_{{dl}}=A{c}_{{dl}}$ is the double layer capacitance of the interface ($F$). Since current generated through faradaic and non-faradaic processes are additive in Eq. 9, it represents parallel paths across the interface, akin to a Randles circuit but with an extra nonlinear charge transfer term that scales with the charge transfer symmetry parameter $\frac{\varepsilon f}{{R}_{{ct}}}.$

To find ${\widetilde{Z}}_{1}\left(\omega \right)$ and ${\widetilde{Z}}_{2}\left(\omega \right)$ for this weakly nonlinear Randles circuit, we insert Eqs. 5 and 6 into Eq. 9, separate the equations into first harmonic (linear) terms multiplied by $\exp \left(j\omega t\right),$ second harmonic terms multiplied by $\exp \left(j2\omega t\right),$ and ignore any higher harmonic terms. With algebraic manipulation, the governing equation for the complex interfacial voltage fluctuations at the excitation frequency $\omega$ is,

$$\begin{eqnarray}&&{(1+{\omega }^{* }j)\widetilde{E}}_{1}={R}_{{ct}}\unicode{x02206}I\end{eqnarray} \tag{ 10 }$$

and the second harmonic interfacial voltage at a frequency of $2\omega$ is

$$\begin{eqnarray}&&{(1+2{\omega }^{* }j)\widetilde{E}}_{2}=-\varepsilon f{\widetilde{E}}_{1}^{2}\end{eqnarray} \tag{ 11 }$$

where ${\omega }^{* }=\omega {R}_{{ct}}{C}_{{dl}}$ is the dimensionless frequency commonly used in Randles circuit analysis. Note that Eq. 11 shows that the second harmonic potential response (${\widetilde{E}}_{2}$) depends on the solution at the excitation frequency (${\widetilde{E}}_{1}$). Algebraic manipulation of Eq. 10, along with the use of Eq. 7, provides the classic linear Randles EIS frequency response at the first harmonic,

$$\begin{eqnarray}&&{\widetilde{Z}}_{1}=\frac{{R}_{{CT}}}{1+{\omega }^{* }j}\end{eqnarray} \tag{ 12a }$$

with real and imaginary parts given as

$$\begin{eqnarray}&&\mathrm{Re}\left({\widetilde{Z}}_{1}\right)={\widetilde{Z}}_{1}^{{\prime} }=\frac{{R}_{{ct}}}{1+{{\omega }^{* }}^{2}}\end{eqnarray} \tag{ 12b }$$

and

$$\begin{eqnarray}&&{Im}({\widetilde{Z}}_{1}){\,=\,\widetilde{Z}}_{1}^{{\rm{\text{'}\text{'}}}}=-\frac{{{\omega }^{* }R}_{{ct}}}{1+{{\omega }^{* }}^{2}}\end{eqnarray} \tag{ 12c }$$

To find the weakly nonlinear 2nd-NLEIS response, Eq. 10 is rearranged and inserted into Eq. 11, then solved for ${\widetilde{E}}_{2}.$ Using Eq. 8 provides the 2nd-NLEIS for the nonlinear Randles circuit,

$$\begin{eqnarray}&&{\widetilde{Z}}_{2}=\frac{-\varepsilon f{R}_{{ct}}^{2}}{1+4{\omega }^{* }j-5{{\omega }^{* }}^{2}-2{{\omega }^{* }}^{3}j}\end{eqnarray} \tag{ 13a }$$

which can be written as real and imaginary terms

$$\begin{eqnarray}&&\mathrm{Re}\left({\widetilde{Z}}_{2}\right)={\widetilde{Z}}_{2}^{{\prime} }=\frac{-\varepsilon f({\widetilde{Z}}_{1}^{{\prime} 2}-{\widetilde{Z}}_{1}^{{\rm{\text{'}\text{'}}}2}+4{\omega }^{* }{\widetilde{Z}}_{1}^{{\prime} }{\widetilde{Z}}_{1}^{{\rm{\text{'}\text{'}}}})}{1+4{{\omega }^{* }}^{2}}\end{eqnarray} \tag{ 13b }$$

and

$$\begin{eqnarray}&&{Im}\left({\widetilde{Z}}_{2}\right)={\widetilde{Z}}_{2}^{{\rm{\text{'}\text{'}}}}=\frac{\varepsilon f\left[2{\omega }^{* }({\widetilde{Z}}_{1}^{{\prime} 2}-{\widetilde{Z}}_{1}^{{\rm{\text{'}\text{'}}}2})-2{\widetilde{Z}}_{1}^{{\prime} }{\widetilde{Z}}_{1}^{{\rm{\text{'}\text{'}}}}\right]}{1+4{{\omega }^{* }}^{2}}\end{eqnarray} \tag{ 13c }$$

Equation 13 shows the size of the 2nd-NLEIS response is directly proportional to the charge transfer asymmetry parameter ( $\varepsilon$ ), with the same characteristic ${R}_{{ct}}{C}_{{dl}}$ time constant as is found in the linear EIS response.

Figure 2 shows Nyquist plots for the first harmonic (linear) EIS and 2nd-NLEIS response of a weakly nonlinear Randles circuit for a series of different charge transfer asymmetry parameters ( $\varepsilon$ ). As expected, the first harmonic linear EIS response in Fig. 2a is a classic Randles semi-circle of width ${R}_{{ct}}$ and characteristic frequency (top of the semi-circle) of $\omega =\frac{1}{{R}_{{ct}}{C}_{{dl}}}.$ The linear EIS response is not sensitive to the charge transfer asymmetry parameter. In contrast, the 2nd-NLEIS signal in Fig. 2b is characterized by a frequency-dependent spiral toward the origin, starting from a real-axis DC limit of ${\widetilde{Z}}_{2}\left(\omega \to 0\right)=-\varepsilon f{R}_{{ct}}^{2}.$ Each 2nd-NLEIS spiral enters three of the four Nyquist plot quadrants. The 2nd-NLEIS response for a half-cell vanishes at all frequencies when charge transfer is symmetric with $\varepsilon =0$.

Zoom In Zoom Out Reset image size

Circuit drawings enable easy visualization of the structure and elements of a linear EIS formulation. Our weakly nonlinear Randles circuit arose from parallel interfacial current paths determined from a two-term nonlinear representation of Butler-Volmer charge transfer (faradaic path) and linear representation of a constant capacitance Helmholtz double layer (non-faradaic path). We illustrate this physics with the nonlinear Randles equivalent circuit drawing shown in Fig. 3, where the charge transfer resistance has a second parallel dashed line to denote the second-order nonlinearity used in the weakly nonlinear expansion. In the two-term expansion of charge-transfer, we now have two physicochemical parameters describing the linear charge transfer resistance and the charge transfer asymmetry. Had we used higher order terms in Eq. 1 and added higher order harmonics to Eq. 6, we would have added extra dashed lines accordingly. A classic capacitor symbol denotes the conventional (linear) non-faradaic current path.

Zoom In Zoom Out Reset image size

As discussed above, ${E}_{n}$ can be treated as zero when all concentrations are time-invariant and uniform, allowing Eq. 4 to be rewritten with common parameters used in analysis:

$$\begin{eqnarray}&&\frac{{d}^{2}}{{R}_{{pore}}{C}_{{dl}}}\frac{{\partial }^{2}E}{\partial {x}^{2}}=\frac{{dE}}{{dt}}+\frac{E}{{C}_{{dl}}{R}_{{ct}}}+\frac{\varepsilon f{E}^{2}}{{C}_{{dl}}{R}_{{ct}}}\end{eqnarray} \tag{ 14 }$$

where ${R}_{{pore}}=\frac{{\rho }_{i}d}{A}$ is the effective electrolyte resistance across the entire porous electrode (${\rm{\Omega }}$), ${C}_{{dl}}=A{S}_{c}{c}_{{dl}}d$ is the double layer capacitance of the porous interface ($F$), ${R}_{{ct}}=\frac{1}{A{S}_{c}{i}_{0}{fd}}$ is the charge transfer resistance of the porous electrode (${\rm{\Omega }}$). Paasch et al. ²⁵ has shown that the continuous form of Eq. 14 can be represented by a transmission-line equivalent circuit, as shown in Fig. 4, though here we have substituted N weakly nonlinear Randles elements for the classic linear transmission line, as elaborated above. In the Supplementary materials, we demonstrate the equivalency of Eq. 14 and Fig. 4 for the two-term weakly nonlinear porous electrode. ²⁸ Likewise, we expect that replacing linear elements with the appropriate nonlinear elements in any sort of sophisticated transmission line model (TLM), such as the 3D TLM proposed by Moškon et al., ^29,30 mono-rail TLM by Solchenbach et al., ³¹ and multi-rail TLM by Siroma et al., ³² should generate the appropriate 2nd-NLEIS signal.

Zoom In Zoom Out Reset image size

Substituting Eq. 6 into Eq. 14, and separating into first and second harmonic terms, again ignoring any higher-order terms, we find a first harmonic governing equation that is identical to the Paasch model (in the limit of zero solid resistivity),

$$\begin{eqnarray}&&\frac{{d}^{2}}{{C}_{{dl}}{R}_{{pore}}}\frac{{\partial }^{2}{\widetilde{E}}_{1}}{\partial {x}^{2}}=\left(j\omega +\frac{1}{{C}_{{dl}}{R}_{{ct}}}\right){\widetilde{E}}_{1}\end{eqnarray} \tag{ 15a }$$

with boundary conditions dictated by the current excitation conditions imposed at the current collector and electrolyte edge of the porous electrode,

$$\begin{eqnarray}&&\frac{\partial {\widetilde{E}}_{1}}{\partial x}{|}_{x=0}=-\frac{{R}_{{pore}}}{d}\unicode{x02206}I\end{eqnarray} \tag{ 15b }$$

and

$$\begin{eqnarray}&&\frac{\partial {\widetilde{E}}_{1}}{\partial x}{|}_{x=d}=0\end{eqnarray} \tag{ 15c }$$

The second harmonic governing equation at twice the excitation frequency ($2\omega$) is written

$$\begin{eqnarray}&&\frac{{d}^{2}}{{C}_{{dl}}{R}_{{pore}}}\frac{{\partial }^{2}{\widetilde{E}}_{2}}{\partial {x}^{2}}=\left(j2\omega +\frac{1}{{C}_{{dl}}{R}_{{ct}}}\right){\widetilde{E}}_{2}+\frac{\varepsilon f{\widetilde{E}}_{1}^{2}}{{C}_{{dl}}{R}_{{ct}}},\end{eqnarray} \tag{ 16a }$$

where we note the non-homogeneous term involving ${\widetilde{E}}_{1}^{2}$ requires the solution to Eq. 15. The homogeneous boundary conditions for the second harmonic are

$$\begin{eqnarray}&&\frac{\partial {\widetilde{E}}_{2}}{\partial x}{|}_{x=0}=0\end{eqnarray} \tag{ 16b }$$

and

$$\begin{eqnarray}&&\frac{\partial {\widetilde{E}}_{2}}{\partial x}{|}_{x=d}=0\end{eqnarray} \tag{ 16c }$$

where Eqs. 16b and 16c ensure that no second harmonic current exists at the current collector or bulk electrolyte boundary. This is a requirement in our impedance formulation because the excitation current has a single frequency, as in Eq. 5. However, the inhomogeneous term in Eq. 16a will generate internal second harmonic currents within the electrode interfaces when charge transfer is asymmetric, i.e., $\varepsilon \ne 0.$

The solution of Eq. 15 is,

$$\begin{eqnarray}&&{\widetilde{E}}_{1}=\unicode{x02206}I{\cdot }{R}_{{pore}}\frac{\cosh \left({\beta }_{1}\frac{\left(d-x\right)}{d}\right)}{{\beta }_{1}\sinh \left({\beta }_{1}\right)}\end{eqnarray} \tag{ 17a }$$

with the parameter

$$\begin{eqnarray}&&{\beta }_{1}{=\left(j\omega {C}_{{dl}}{R}_{{pore}}+\frac{{R}_{{pore}}}{{R}_{{ct}}}\right)}^{\frac{1}{2}}\end{eqnarray} \tag{ 17b }$$

The linear EIS response for the weakly nonlinear porous electrode is found by taking Eq. 17 at the electrolyte interface ($x=0$) and dividing by the current excitation amplitude, i.e.,

$$\begin{eqnarray}&&{\widetilde{Z}}_{1}\left(\omega \right)=\frac{{R}_{{pore}}\coth \left({\beta }_{1}\right)}{{\beta }_{1}},\end{eqnarray} \tag{ 18 }$$

which is identical to the linear Paasch model in the limit of zero solid resistivity.

Equation 17b shows there are two relevant time constants in this porous electrode formulation, yielding two dimensionless frequencies, ${\omega }^{* }=\omega {C}_{{dl}}{R}_{{ct}}$ and ${\omega }^{\dagger }=\omega {C}_{{dl}}{R}_{{pore}}.$ The ratio between these frequencies ($\frac{{R}_{{pore}}}{{R}_{{ct}}}$) dictates how accessible the electrode pore structure is to current, namely, if it is a 'thin' uniformly accessible electrode (${R}_{{pore}}\ll {R}_{{ct}}$), a "thick" electrode with a poorly accessible interfacial area (${R}_{{pore}}\gg {R}_{{ct}}$), or somewhere in-between. The D.C. limit for the electrode impedance makes the implications of 'thin' or "thick" limits clearer on the width of the impedance arc:

$$\begin{eqnarray}&&{\widetilde{Z}}_{1}\left(\omega \to 0\right)=\frac{{R}_{{pore}}}{\sqrt{\frac{{R}_{{pore}}}{{R}_{{ct}}}}}\coth {\rm{}}\left(\sqrt{\frac{{R}_{{pore}}}{{R}_{{ct}}}}\right)\end{eqnarray} \tag{ 19 }$$

The solution for Eq. 16a is

$$\begin{eqnarray*}&&{\widetilde{E}}_{2}={\widetilde{C}}_{1}\sinh \left({\beta }_{2}\frac{x}{d}\right)+{\widetilde{C}}_{2}\cosh \left({\beta }_{2}\frac{x}{d}\right)\end{eqnarray*}$$

$$\begin{eqnarray}&&-\frac{\varepsilon f{R}_{{pore}}^{3}{\left(\unicode{x02206}I\right)}^{2}}{{R}_{{ct}}\,{{(\beta }_{1}\sinh \left({\beta }_{1}\right))}^{2}}\left[\frac{\cosh \left(2{\beta }_{1}\frac{\left(d-x\right)}{d}\right)}{2\left({\beta }_{2}-2{\beta }_{1}\right)\left({\beta }_{2}+2{\beta }_{1}\right)}+\frac{1}{2{\beta }_{2}^{2}}\right]\end{eqnarray} \tag{ 20a }$$

where

$$\begin{eqnarray}&&{\widetilde{C}}_{1}=-\frac{\varepsilon f{R}_{{pore}}^{3}{\left(\unicode{x02206}I\right)}^{2}}{{R}_{{ct}}\,{{(\beta }_{1}\sinh \left({\beta }_{1}\right))}^{2}}\begin{array}{c}\left[\frac{{\beta }_{1}\sinh \left(2{\beta }_{1}\right)}{{\beta }_{2}\left({\beta }_{2}-2{\beta }_{1}\right)\left({\beta }_{2}+2{\beta }_{1}\right)}\right],\end{array}\end{eqnarray} \tag{ 20b }$$

$$\begin{eqnarray}&&{\widetilde{C}}_{2}=-{\widetilde{C}}_{1}\coth \left({\beta }_{2}\right),\end{eqnarray} \tag{ 20c }$$

and

$$\begin{eqnarray}&&{\beta }_{2}={\left(j2\omega {C}_{{dl}}{R}_{{pore}}+\frac{{R}_{{pore}}}{{R}_{{ct}}}\right)}^{\frac{1}{2}}\end{eqnarray} \tag{ 20d }$$

Setting Eq. 20a to its value at $x=0,$ the electrolyte boundary, and dividing by ${\left(\unicode{x02206}I\right)}^{2}$ yields the 2nd-NLEIS response:

$$\begin{eqnarray*}&&{\widetilde{Z}}_{2}\left(\omega \right)=\frac{\varepsilon f{R}_{{pore}}^{3}}{{R}_{{ct}}{{(\beta }_{1}\sinh \left({\beta }_{1}\right))}^{2}}\end{eqnarray*}$$

$$\begin{eqnarray}&&\begin{array}{c}\left[\left[\frac{{\beta }_{1}\sinh \left(2{\beta }_{1}\right)}{{\beta }_{2}\left({\beta }_{2}-2{\beta }_{1}\right)\left({\beta }_{2}+2{\beta }_{1}\right)}\right]\coth \left({\beta }_{2}\right)\right.\\ \left.-\,\left[\frac{\cosh \left(2{\beta }_{1}\right)}{2\left({\beta }_{2}-2{\beta }_{1}\right)\left({\beta }_{2}+2{\beta }_{1}\right)}+\frac{1}{2{\beta }_{2}^{2}}\right]\right]\end{array}.\end{eqnarray} \tag{ 21 }$$

The functional forms of ${\beta }_{1}$ and ${\beta }_{2}$ suggest the 2nd-NLEIS response shares the same dimensionless frequencies as the linear system, namely, ${\omega }^{* }=\omega {C}_{{dl}}{R}_{{ct}}$ and ${\omega }^{\dagger }=\omega {C}_{{dl}}{R}_{{pore}}.$ The magnitude for the 2nd-NLEIS spiral is set by the D.C. limit,

$$\begin{eqnarray*}&&{\widetilde{Z}}_{2}\left(\omega \to 0\right)\,=\end{eqnarray*}$$

$$\begin{eqnarray}&&-\,\frac{\varepsilon f{R}_{{ct}}{R}_{{pore}}}{{\sinh }^{2}\left(\sqrt{\frac{{R}_{{pore}}}{{R}_{{ct}}}}\right)\,}\left[1+\frac{1}{3}{\sinh }^{2}\left(\sqrt{\frac{{R}_{{pore}}}{{R}_{{ct}}}}\right)\right].\end{eqnarray} \tag{ 22 }$$

Figure 5 shows the Nyquist plots for the first harmonic (linear) EIS response given by Eq. 18 and the 2nd-NLEIS response given by Eq.  21 for a "thick" porous electrode where ${R}_{{pore}}=10{R}_{{ct}}.$ The first harmonic response shown in Fig. 5a is characterized by a 'semi-teardrop' shape that is characteristic for thick electrodes with poor utilization of internal interfacial area. In the thick electrode limit, the frequency at the top of the EIS curve (maximum imaginary impedance) is ${\omega }^{* }=\sqrt{3}.$ Though not shown, 'thin' porous electrodes (${R}_{{pore}}\ll {R}_{{ct}}$) retain the classic shape of a Randles semi-circle except at extremely high frequencies. The 2nd-NLEIS response shown in Fig. 5b has a spiral-like character akin to the weakly nonlinear Randles circuit, but with subtle differences in shape. An example of the subtle shape difference is that the 2nd-NLEIS curves in Fig. 5b only traverse two of the four Nyquist plot quadrants, whereas the Randles case in Fig. 2b had a three-quadrant spiral. As is always true for the 2nd-NLEIS, the half-cell response is proportional to the charge transfer asymmetry parameter and has a null response for symmetric charge transfer.

Zoom In Zoom Out Reset image size

Diffusion dynamics in the solid or electrolyte phase of an electrochemical system can give rise to a weakly nonlinear Warburg-like diffusion impedance. Here we expand the analysis above to include solid-state bounded diffusion in the nonlinear Randles and porous electrode impedance responses. The archetypical problem considered here involves Li-ion insertion materials, where the equilibrium potential is a function of the concentration of lithium ions inserted into the solid electrode (denoted ${c}_{s}$). A Taylor series expansion can be used to define the second order processes that describe the nonlinear oscillating equilibrium potential (${\widetilde{E}}_{n}$) of a fully relaxed non-hysteretic system,

$$\begin{eqnarray}&&{\widetilde{E}}_{n}={\left.\left(\frac{\partial {E}_{n}}{\partial {c}_{s}}\right)\right|}_{{SoC}}{\widetilde{c}}_{s}+\,\frac{1}{2}{\left.\left(\frac{{\partial }^{2}{E}_{n}}{\partial {c}_{s}^{2}}\right)\right|}_{{SoC}}{\widetilde{c}}_{s}^{2}+h.o.t\end{eqnarray} \tag{ 23 }$$

where ${\widetilde{c}}_{s}$ is the oscillating concentration of inserted ions at the solid interface, ${\left.\left(\frac{\partial {E}_{n}}{\partial {c}_{s}}\right)\right|}_{{SoC}}$ and ${\left.\left(\frac{{\partial }^{2}{E}_{n}}{\partial {c}_{s}^{2}}\right)\right|}_{{SoC}}$ are the slope and curvature of the equilibrium potential at the given state of charge (SoC) for the insertion material. To be consistent with prior results and the two-term expansion in Eq. 23, the weakly nonlinear solid state surface concentration fluctuations are expected to follow

$$\begin{eqnarray}&&{\widetilde{c}}_{s}={\widetilde{c}}_{s,1}\exp \left(j\omega t\right)+{\widetilde{c}}_{s,2}\exp \left(j2\omega t\right)\end{eqnarray} \tag{ 24 }$$

where ${\widetilde{c}}_{s,1}$ and ${\widetilde{c}}_{s,2}$ are the first and second harmonic coefficients for the fluctuating insertion concentrations at the solid surface. Substituting Eqs. 24 into 23 and ignoring higher harmonics yields the weakly nonlinear oscillating equilibrium potential:

$$\begin{eqnarray*}&&{\widetilde{E}}_{n}={\left.\left(\frac{\partial {E}_{n}}{\partial {c}_{s}}\right)\right|}_{{SoC}}{\widetilde{c}}_{s,1}\exp \left(j\omega t\right)\end{eqnarray*}$$

$$\begin{eqnarray}&&+\,\left[{\left.\left(\frac{\partial {E}_{n}}{\partial {c}_{s}}\right)\right|}_{{SoC}}{\widetilde{c}}_{s,2}+\frac{1}{2}{\left.\left(\frac{{\partial }^{2}{E}_{n}}{\partial {c}_{s}^{2}}\right)\right|}_{{SoC}}{\widetilde{c}}_{s,1}^{2}\right]\exp \left(j2\omega t\right)\end{eqnarray} \tag{ 25 }$$

An oscillating faradaic current of the form

$$\begin{eqnarray}&&{I}_{f}={\widetilde{I}}_{f,1}\exp \left(j\omega t\right)+{\widetilde{I}}_{f,2}\exp \left(j2\omega t\right)\end{eqnarray} \tag{ 26 }$$

is consistent with our prior solutions for both planar and porous electrodes, even though the total current represented by Eq. 5 has no second harmonic component entering or leaving the electrode. The transfer function relating oscillating surface concentration to oscillating faradaic current is known as the diffusion impedance (${\widetilde{Z}}_{D}$), which we write as ³³

$$\begin{eqnarray}&&{\widetilde{Z}}_{D,1}\left(\omega \right)={\widetilde{Z}}_{D}(\omega )={\left.\left(\frac{\partial {E}_{n}}{\partial {c}_{s}}\right)\right|}_{{SoC}}\frac{{\left.{\widetilde{c}}_{s,1}\right|}_{x=l}}{{\widetilde{I}}_{f,1}}\end{eqnarray} \tag{ 27a }$$

and

$$\begin{eqnarray}&&{\widetilde{Z}}_{D,2}\left(\omega \right)={\widetilde{Z}}_{D}\left(2\omega \right)={\left.\left(\frac{\partial {E}_{n}}{\partial {c}_{s}}\right)\right|}_{{SoC}}\frac{{\left.{\widetilde{c}}_{s,2}\right|}_{x=l}}{{\widetilde{I}}_{f,2}}.\end{eqnarray} \tag{ 27b }$$

Inserting Eq. 27 into Eq. 25 produces a compact form of the oscillating equilibrium potential

$$\begin{eqnarray*}&&{\widetilde{E}}_{n}={\widetilde{Z}}_{D,1}{\widetilde{I}}_{f,1}\exp \left(j\omega t\right)\end{eqnarray*}$$

$$\begin{eqnarray}&&+\,\left[{\widetilde{Z}}_{D,2}{\widetilde{I}}_{f,2}+\kappa {\left({\widetilde{Z}}_{D,1}{\widetilde{I}}_{f,1}\right)}^{2}\right]\exp \left(j2\omega t\right)\end{eqnarray} \tag{ 28 }$$

where we have defined a new differential voltage parameter $\kappa$ (${V}^{-1}$)

$$\begin{eqnarray}&&\kappa =\frac{{\left.\left(\frac{{\partial }^{2}{E}_{n}}{\partial {c}_{s}^{2}}\right)\right|}_{{SoC}}}{2{\left.{\left(\frac{\partial {E}_{n}}{\partial {c}_{s}}\right)}^{2}\right|}_{{SoC}}}\end{eqnarray} \tag{ 29 }$$

that scales with the curvature of the open circuit potential at a given SoC. Charge transfer-driven diffusion nonlinearities occur at frequencies relevant to diffusive timescales set via ${\widetilde{Z}}_{D}.$ Expressions for the bounded diffusion impedance ${\widetilde{Z}}_{D}(\omega )$ have been calculated for a variety of insertion material geometries and are summarized in Table I. ^34–36

Table I. Bounded diffusion impedance for different insertion material geometries. The diffusion timescale $\tau =\tfrac{{l}^{2}}{D}$ ($s$) is tied to the insertion material characteristic dimension $l(m)$ and inserted ion diffusion coefficient $D({m}^{2}/s)$. Bounded diffusion impedance scales with ${A}_{w}=\left(-\tfrac{\partial {E}_{n}}{\partial {c}_{s}}\right)\left(\tfrac{l}{nFD{A}_{i}}\right),$ where ${A}_{i}$ is the interfacial area of the electrode. The functions ${I}_{0}$ and ${I}_{1}$ denote the modified Bessel functions of zero and first order, respectively.

	Planar	Cylindrical	Spherical
${\widetilde{Z}}_{D}\left(\omega \right)$	$\frac{{A}_{w}\coth \left(\sqrt{j\omega \tau }\right)}{\sqrt{j\omega \tau }}$	$\frac{{A}_{w}{I}_{0}(\sqrt{j\omega \tau })}{\sqrt{j\omega \tau }{I}_{1}(\sqrt{j\omega \tau })}$	$\frac{{A}_{w}\tanh \left(\sqrt{j\omega \tau }\right)}{\sqrt{j\omega \tau }-\tanh \left(\sqrt{j\omega \tau }\right)}$

The first and second harmonic faradaic current modulations can be found through algebraic manipulation of Eqs. 1, 6, 26, and 28 to identify the coefficients multiplying $\exp \left(j\omega t\right),$ that is,

$$\begin{eqnarray}&&{I}_{f,1}=\frac{{\widetilde{E}}_{1}}{{\widetilde{Z}}_{D,1}+{R}_{{ct}}}\end{eqnarray} \tag{ 30a }$$

and, likewise, the coefficient multiplying $\exp \left(j2\omega t\right)$ is:

$$\begin{eqnarray}&&{I}_{f,2}=\frac{\left[{\widetilde{E}}_{2}-\kappa {\left({\widetilde{Z}}_{D,1}{I}_{f,1}\right)}^{2}\right]+\varepsilon f{\left[{\widetilde{E}}_{1}-{\widetilde{Z}}_{D,1}{I}_{f,1}\right]}^{2}}{{\widetilde{Z}}_{D,2}+{R}_{{ct}}}\end{eqnarray} \tag{ 30b }$$

These expressions hold for both the planar and the porous electrode.

Equation 3 provides the template for treating a weakly nonlinear planar electrode when faradaic charge transfer is accompanied by diffusion of inserted ions, that is, when Eqs. 28 and 30 are used with Eq. 3. The resulting nonlinear equivalent circuit diagram is represented by Fig. 6.

Zoom In Zoom Out Reset image size

Substituting Eqs. 6, 28 , and 30 into Eq. 3 followed by algebraic manipulation and separation of equations into first and second harmonics produces expressions for ${\widetilde{E}}_{1}$ and ${\widetilde{E}}_{2}.$ Equation 7 is used with the ${\widetilde{E}}_{1}$ solution to write the first harmonic (linear) EIS response

$$\begin{eqnarray}&&{\widetilde{Z}}_{1}=\frac{{R}_{{ct}}}{\frac{{R}_{{ct}}}{{R}_{{ct}}+{\widetilde{Z}}_{D,1}}+j{\omega }^{* }}\end{eqnarray} \tag{ 31 }$$

while solutions to ${\widetilde{E}}_{2}$ are used with Eq. 8 to write the second harmonic NLEIS response

$$\begin{eqnarray}&&{\widetilde{Z}}_{2}=\frac{{R}_{{ct}}}{\left(j2{\omega }^{* }+\frac{{R}_{{ct}}}{{\widetilde{Z}}_{D,2}+{R}_{{ct}}}\right)}\frac{\left[\kappa {\left[\frac{{\widetilde{Z}}_{D,1}}{{\widetilde{Z}}_{D,1}+{R}_{{ct}}}\right]}^{2}-\varepsilon f{\left[\frac{{R}_{{ct}}}{{\widetilde{Z}}_{D,1}+{R}_{{ct}}}\right]}^{2}\right]}{{\widetilde{Z}}_{D,2}+{R}_{{ct}}}{\left(\frac{{R}_{{ct}}}{\frac{{R}_{{ct}}}{{R}_{{ct}}+{\widetilde{Z}}_{D,1}}+j{\omega }^{* }}\right)}^{2}\end{eqnarray} \tag{ 32 }$$

Equations 31 and 32 enable the simulation of the full spectrum linear EIS and 2nd-NLEIS response for a planar thin film of insertion material on a current collector substrate. Figure 7a shows a Nyquist plot for linear EIS, where a Randles charge transfer semi-circle at higher frequencies is seen to merge into a 45-degree diffusive impedance at intermediate frequencies and, ultimately, the chemical capacitance of the bounded thin film is observed at sufficiently low frequencies. The corresponding 2nd-NLEIS response is explored as a function of the charge transfer asymmetry parameter ($\varepsilon$) in Fig. 7b and as a function of the voltage curvature parameter ($\kappa$) in Fig. 7c. Figure 7b shows the same three-quadrant 2nd-NLEIS charge transfer spiral as the weakly nonlinear Randles circuit presented in Fig. 2, but as the frequency becomes low enough to activate diffusive processes, the weakly nonlinear bounded diffusion impedance terms in Eq. 32 becomes apparent. Figure 7c shows the effect of $\kappa$ on the nonlinear second harmonic response. As shown in the Nyquist plot, the direction and magnitude of the weakly nonlinear diffusion response are governed by the $\kappa .$ For the parameters shown here, the case of $\kappa$ equal to zero is close to the pure charge transfer behavior with Eq. 32 reverting to the pure charge transfer DC limit. One could imagine when the charge transfer time scale and diffusion time scale are comparable to each other, a larger deviation from pure charge transfer reaction might displayed even when $\kappa$ is zero. The direction-sensitive features of 2nd-NLEIS — where realistic charge transfer and diffusion parameters can produce signals in all four quadrants of a Nyquist plot — significantly boost the mathematical discriminating power of 2nd-NLEIS over conventional one-quadrant linear EIS. We discuss this parameter sensitivity as an especially important feature for full-cell studies in the Implications and Concluding Remarks section.

Zoom In Zoom Out Reset image size

Equation 4 is the starting point for describing a porous insertion electrode with solid-state diffusive transport. The corresponding weakly nonlinear equivalent circuit representation of such a porous electrode is shown in Fig. 8.

Zoom In Zoom Out Reset image size

Substituting Eqs. 5, 28 , and 30 into Eq. 4 and separating terms into first and second harmonic equations, while ignoring higher order harmonics, yields the first harmonic (linear) governing equation:

$$\begin{eqnarray}&&\frac{{d}^{2}}{{C}_{{dl}}{R}_{{pore}}}\frac{{\partial }^{2}{\widetilde{E}}_{1}}{\partial {x}^{2}}=\left[j\omega +\frac{1}{{C}_{{dl}}\left({\widetilde{Z}}_{D,1}+{R}_{{ct}}\right)}\right]{\widetilde{E}}_{1}\end{eqnarray} \tag{ 33 }$$

with identical boundary conditions to Eq. 15b and 15c. The solution to Eq. 33 is identical to Eq. 17a, but with Eq. 17b replaced by the diffusion impedance modified parameter

$$\begin{eqnarray}&&{\beta }_{1}^{D}{=\left(j\omega {C}_{{dl}}{R}_{{pore}}+\frac{{R}_{{pore}}}{{\widetilde{Z}}_{D,1}+{R}_{{ct}}}\right)}^{\frac{1}{2}}\end{eqnarray} \tag{ 34 }$$

yielding the linear EIS impedance

$$\begin{eqnarray}&&{\widetilde{Z}}_{1}^{D}\left(\omega \right)=\frac{{R}_{{pore}}\coth \left({\beta }_{1}^{D}\right)}{{\beta }_{1}^{D}}\end{eqnarray} \tag{ 35 }$$

Formulating the second harmonic porous electrode response when solid state bounded diffusion effects are included yields a governing equation and boundary conditions identical in form to Eq. 16a, resulting in a solution having an identical form to Eq. 20a and 2nd-NLEIS response

$$\begin{eqnarray*}&&{\widetilde{Z}}_{2}^{D}\left(\omega \right)=\frac{{R}_{{pore}}^{3}{k}^{D}}{{R}_{{ct}}{\left[{\beta }_{1}^{D}\sinh \left({\beta }_{1}^{D}\right)\right]}^{2}}\end{eqnarray*}$$

$$\begin{eqnarray}&&\begin{array}{c}\left[\left[\frac{{\beta }_{1}^{{\rm{D}}}\sinh \left(2{\beta }_{1}^{D}\right)}{{\beta }_{2}^{D}\left({\beta }_{2}^{D}-2{\beta }_{1}^{D}\right)\left({\beta }_{2}^{D}+2{\beta }_{1}^{D}\right)}\right]\coth \left(\Space{0ex}{2.5ex}{0ex},(,{\beta }_{2}^{D},\Space{0ex}{2.5ex}{0ex}\right)\right.\\ \left.-\,\left[\frac{\cosh \left(2{\beta }_{1}^{D}\right)}{2\left({\beta }_{2}^{D}-2{\beta }_{1}^{D}\right)\left({\beta }_{2}^{D}+2{\beta }_{1}^{D}\right)}+\frac{1}{2{{\beta }_{2}^{D}}^{2}}\right]\right]\end{array}\end{eqnarray} \tag{ 36 }$$

with the diffusion impedance modified parameters

$$\begin{eqnarray}&&{\beta }_{2}^{D}{=\,\left(j2\omega {C}_{{dl}}{R}_{{pore}}+\frac{{R}_{{pore}}}{{\widetilde{Z}}_{D,2}+{R}_{{ct}}}\right)}^{\frac{1}{2}}\end{eqnarray} \tag{ 37a }$$

and

$$\begin{eqnarray}&&{k}^{D}=\frac{\left[-{R}_{{ct}}\kappa {\left[\frac{{\widetilde{Z}}_{D,1}}{{\widetilde{Z}}_{D,1}+{R}_{{ct}}}\right]}^{2}+{R}_{{ct}}\varepsilon f{\left[\frac{{R}_{{ct}}}{{\widetilde{Z}}_{D,1}+{R}_{{ct}}}\right]}^{2}\right]}{{\widetilde{Z}}_{D,2}+{R}_{{ct}}}.\end{eqnarray} \tag{ 37b }$$

being the only differentiating factor between Eqs.  36 and  21.

Equations 35 and  36 enable the simulation of the full spectrum linear EIS and 2nd-NLEIS response for a porous insertion electrode on a current collector substrate. Figure 9a shows a Nyquist plot for linear EIS, where a typical "thick" porous electrode charge transfer teardrop shape is seen at frequencies above 100 mHz, akin to what is seen in Fig. 5a without any diffusion effects. At frequencies below 100 mHz, the effect of diffusive transport becomes evident as the 45-degree diffusive impedance transitions to a chemical capacitance behavior in the bounded thin film at the lowest frequencies shown. Here we assume the porous electrode was built from planar (platelet-like) particles, hence the "planar" functional form is used for the diffusion impedance in Table I. The corresponding 2nd-NLEIS response is explored as a function of the charge transfer asymmetry parameter ($\varepsilon$) in Fig. 9b and as a function of the potential curvature parameter ($\kappa$) in Fig. 9c. Figure 9b shows the same two-quadrant 2nd-NLEIS charge transfer spiral at frequencies above 100 mHz as we saw in Fig. 5, but as the frequency becomes lower, the weakly nonlinear bounded diffusion impedance terms in Eqs.  36 and 37 becomes apparent. Figure 9c shows the effect of $\kappa$ on the scaling and quadrants where the low frequency nonlinear second harmonic diffusion processes appear. In contrast to the planar electrode, there is a distinct weakly nonlinear diffusion response for the porous electrode even when $\kappa$ equals to zero, which is primarily due to the inhomogeneity of the current distribution with the presence of pore electrolyte.

Zoom In Zoom Out Reset image size

Linearization used in traditional EIS techniques leads to half-cell model degeneracy with both physics-based models ¹² or equivalent circuit models. ¹¹ The challenge of half-cell degeneracy is compounded when analyzing two-electrode cells. The origin of this compounding impact is that a two-electrode cell has an EIS response that is the sum of half-cell EIS signals; thus, two-electrode linear EIS is a "positive parity" signal owing to the additive trait of half-cells. This positive parity trait makes it even more difficult to assign the physicochemical contributions to the EIS signal of a particular electrode.

2nd-NLEIS is effective at addressing the challenge of model degeneracy in half-cell measurements. ¹² Perhaps more importantly, the two-electrode 2nd-NLEIS signal arises from the difference between the half-cell 2nd-NLEIS responses; ²² thus, two electrode 2nd-NLEIS is a "negative parity" signal owing to subtractive traits of half-cells. The complementary parity between EIS and 2nd-NLEIS, when analyzed with a common physics-based model, makes assigning half-cell physicochemical processes from two-electrode measurements much more feasible.

To demonstrate the importance of using complementary parity signals to assign physicochemical process, we use the half-cell modeling above to build four different two-electrode models that each fairly closely represent the linear EIS response from two-electrode battery data presented by Murbach et al. ²² The four distinct models were built by simply "swapping" charge transfer and diffusion-related physicochemical parameters systematically among the positive and negative electrodes, creating the four parameter sets summarized in Table II. Figure 10a shows the simulated Nyquist plots for each two-electrode linear EIS model. The positive parity additive response makes it difficult to assign the charge transfer and diffusion properties to one electrode or another. Though not truly degenerate (there are subtle differences among the four cases), the two-electrode EIS responses are practically indistinguishable. The Nyquist plots for the four corresponding two-electrode 2nd-NLEIS models are shown in Fig. 10b. In contrast to Fig. 10a, each negative-parity NLEIS model for the two-electrode cells has a highly distinct response that is sensitive to the assignment of specific physicochemical parameters to specific electrodes. In short, the complementary parity between linear EIS and 2nd-NLEIS in a two-electrode configuration provides the discriminating power needed to assign half-cell processes in a two-electrode configuration that lacks a separate reference measurement.

Table II. Parameters of first and second harmonic impedance simulation for four nearly degenerated linear first harmonic models and their corresponding nonlinear second harmonic models.

Physicochemical Parameter			Model 1	Model 2	Model 3	Model 4
Cathode	Charge Transfer	${R}_{{pore},c}\,[{\rm{\Omega }}]$	$0.0198$	$0.0198$	$0.00860$	$0.00860$
		${R}_{c}\,[{\rm{\Omega }}]$	$0.00749$	$0.00749$	$0.00214$	$0.00214$
		${C}_{c}\,[F]$	$5.03$	$5.03$	$0.390$	$0.390$
	Diffusion	${A}_{W,c}\,[{\rm{\Omega }}]$	$0.138$	$0.143$	$0.138$	$0.143$
		${\tau }_{c}[s]$	$3547$	$13595$	$3547$	$13595$
	NLEIS Parameters	${\kappa }_{c}\,[{V}^{-1}]$	$9.05$	$9.05$	$39.5$	$39.5$
		${\varepsilon }_{c}[]$	$0.0598$	$0.0598$	$0$	$0$
Anode	Charge Transfer	${R}_{{pore},a}\,[{\rm{\Omega }}]$	$0.00860$	$0.00860$	$0.0198$	$0.0198$
		${R}_{a}\,[{\rm{\Omega }}]$	$0.00214$	$0.00214$	$0.00749$	$0.00749$
		${C}_{a}\,[F]$	$0.390$	$0.390$	$5.03$	$5.03$
	Diffusion	${A}_{W,a}\,[{\rm{\Omega }}]$	$0.143$	$0.138$	$0.143$	$0.138$
		${\tau }_{a}[s]$	$13595$	$3547$	$13595$	$3547$
	NLEIS Parameters	${\kappa }_{a}\,[{V}^{-1}]$	$39.5$	$39.5$	$9.05$	$9.05$
		${\varepsilon }_{a}[]$	$0$	$0$	$0.0598$	$0.0598$

Zoom In Zoom Out Reset image size

The implications for the half-cell and two-electrode cell analysis given here, when combined with our prior work on the experimental accessibility of two-harmonic measurements, is that 2nd-NLEIS is a truly complementary measurement that should be added to the EIS toolbox in order to take advantage of the discriminating power described above. To build confidence in 2nd-NLEIS based methods, adoption of the advanced quality assurance and analysis methods that have become fairly routine in the linear EIS community ⁹ need to be advanced to the second harmonic. A key question that underpins the efficient adoption of 2nd-NLEIS is how to systematically test that an experimental modulation is adequate to excite the second-order weakly nonlinear effects analyzed here, without over-driving the system into third-order (and higher) nonlinear processes that will distort the linear EIS signals, avoiding the need to use more sophisticated series analysis to extract any meaningful parameters (e.g., see analysis in Murbach et al.'s work ²³ ). In the experiments we've reported previously, modulated voltages of roughly 10–15 mV (at frequencies that activate charge transfer and diffusion processes) appear to nicely drive second-order (two term, two harmonic) responses without overly exciting higher-order distortions. ²²

As shown above, building the second harmonic nonlinear impedance response from two-term (weakly nonlinear) expansions of interfacial and transport models creates a mathematical structure that is consistent with the parallel and series circuit element additivity rules of linear circuits. We formally introduced the consequences of two-term expansions of the Butler-Volmer charge transfer model and solid-state bounded diffusion. In contrast, we chose a constant coefficient Helmholtz model for the double layer which avoids nonlinear contributions to the non-faradaic current; this allowed us to show that mixing linear and nonlinear elements is easily accomplished and it also aligned the work here with previous nonlinear pseudo-2-D modeling of porous insertion electrodes. ²³

For completeness, we note that more general descriptions of the electrochemical double layer require a replacement of Eq. 2, given that interfacial capacitance is no longer a constant:

$$\begin{eqnarray}&&{{\rm{i}}}_{{nf}}\,=\,{{\rm{c}}}_{{dl}}\frac{{dE}}{{dt}}+E\frac{d{{\rm{c}}}_{{dl}}}{{dt}}\,\end{eqnarray} \tag{ 38a }$$

Most double-layer models represent the double layer capacitance as a function of interfacial potential $E.$ Because 2nd-NLEIS perturbs the interface from its equilibrium state with small but finite modulations, it is reasonable to use a two-term Taylor series expansion for the capacitance, yielding a universal "local" form for the capacitance of

$$\begin{eqnarray}&&{c}_{{dl}}={c}_{{dl},n}\left[1+{\rm{\Theta }}f(E-{E}_{n})\right],\end{eqnarray} \tag{ 38b }$$

where ${c}_{{dl},n}$ is the capacitance of the interface at the equilibrium potential and ${\rm{\Theta }}$ is a model-dependent nondimensional constant. The combination of Eqs. 38a and 38b creates a direct replacement for Eq. 2 that introduces second order nonlinearity into the governing equation for non-faradaic currents. Had Eq. 38 been used throughout this paper, the linear EIS response would be unchanged, but the 2nd-NLEIS response would have many extra terms and an extra parameter ${\rm{\Theta }}$ in the solutions for each specific case studied.

The selection of a particular double layer structure yields a detailed interpretation of the model-dependent constant ${\rm{\Theta }}.$ For example, a Gouy-Chapman double layer with monovalent ions yields ${\rm{\Theta }}=0.5\,{\rm{* }}\,\tanh \left[0.5f\left({E}_{n}-{E}_{{pzc}}\right)\right],$ where ${E}_{{pzc}}$ is the potential of zero charge (PZC) for the interface. For the Gouy-Chapman model, interfacial capacitance has a local minimum at the PZC with ${\rm{\Theta }}=0,$ resulting in no 2nd-NLEIS contribution from the double layer. More broadly, nonlinear double layer contributions to 2nd-NLEIS spectra are likely to only be important at highly polarizable interfaces and at high frequencies where faradaic charge transfer contributions are negligible. In short, though not used here, non-faradaic currents can contribute to 2nd-NLEIS through the second-order nonlinearity generated by a non-constant double layer capacitance. The form of the contribution, Eq. 38b, is universal but the interpretation of the new nonlinear parameter ${\rm{\Theta }}$ is dependent on the interfacial structure of models. For electrochemical energy systems we mostly study, facile charge transfer is a key feature of the electrochemical interface, so the impact of weak double-layer nonlinearity is diminished.

Shown in Table III is a visual representation of the linear and two-term weakly nonlinear elements that underpin all the analysis above. We adapt traditional linear EIS circuit element representations and add a parallel dashed line to represent the second-order nonlinear form of each linear circuit element. Adding a parallel dashed line conveys the use of a two-term expansion to represent the second-order nonlinear electrochemical processes being analyzed. As shown in Figs. 3, 4, 6, and 8, one can mix and match linear and second-order impedance elements to visualize the mix of physical assumptions that underpin the linear and second harmonic response of the corresponding Nyquist plots. Such development of visual representation emphasizes the ability to substitute linear circuit elements with appropriate nonlinear element circuits to take advantage of the existing TLM models ^{29–32,37,38} in describing different electrochemical processes, while obtaining additional insights about reaction kinetics and thermodynamics with the complementary 2nd-NLEIS analysis.

Table III. Circuit element representations for the first and second harmonic impedance response when a two-term expansion is used to represent a nonlinear electrochemical process.

Circuit Elements	Linear EIS	2nd Harmonic NLEIS
Resistor
Capacitor
Diffusion element

To conclude, second harmonic nonlinear electrochemical impedance spectroscopy (2nd-NLEIS) is a powerful tool for probing complementary information to traditional electrochemical impedance spectroscopy (EIS), extracting meaningful non-linear contributions to charge transfer, mass transport, and thermodynamics of the electrochemical systems. However, the application of the 2nd-NLEIS has been hampered by the lack of simple models able to extract meaningful parameters. To resolve this problem, we systematically developed simple physics-based models for planar electrochemical interfaces and porous electrodes that are governed by a Helmholtz double layer, Butler–Volmer kinetics, and solid-state bounded Fickian diffusion. The limiting cases and dependences of the second harmonic impedance response on charge transfer asymmetry and thermodynamics are explored. The importance of using the complementary parity between linear EIS and 2nd-NLEIS as a strategy to separate half-cell parameters for two-electrode measurements is discussed. For completeness, the nonlinear double-layer response for non-faradaic current is briefly introduced for the interested reader. Lastly, we try to depict a sensible visual representation of second order nonlinear circuit elements as a shorthand way to identify the mathematical underpinning and structure of the analysis, helping bridge the gap between model-based analysis, experiments, and the intuition needed for the physical interpretation of EIS and 2nd-NLEIS data.