Lithium-ion Batteries Solid-phase Diffusion Approximations, Transient Response, Electrolyte Concentrations delayed in Time Stability Analysis under Parameters Variation

Lithium – ion batteries are described by mathematical model which involve variables such as electrolyte concentrations which are volume averaged over respective region. Other possible dependent variables can be solution phase potential and solution phase concentration. The schematic Lithium – ion cell sandwich defines three important distances $l_p$ , $l_s$ , and $l_n$ for positive electrode length, separator length, and negative electrode length, respectively. The perturbation techniques define the model and the number of solved equations. Lithium – ion batteries are characterized by three regions with differ physical properties and characterized equations in each region. The Lithium – ion model characterized ordinary differential equations (ODEs) are converted to delay differential equations (DDEs) which reflect the electrolyte concentrations delayed in time due to inconsistency and interferences in the solid – phase diffusion approximations and transient response. The key model delayed variables are electrolyte concentrations which are volume averaged over respective region

( $C_{\text{ave}}^{\text{cathode}}$ , $C_{\text{ave}}$ , $C_{\text{ave}}^{\text{anode}}$ ). The three DDEs Lithium – ion model gives the prediction for electrolyte concentrations under influence of variables delay in time. The standard local stability analysis about any one of the equilibrium points of Lithium – ion batteries model is inspected for stability and stability switching behaviour under delay parameters in time variation.

Export Citation CrossMark Publons Harvard Library HOLLIS GrowKudos Search IT Google Scholar Academic Microsoft Scilit Semantic Scholar Universite de Paris UW Libraries SJSU King Library NUS Library McGill DET KGL BIBLiOTEK JCU Discovery Universidad De Lima WorldCat DTU VU on WorldCat ResearchGate

A lithium – ion battery is unit which store energy and is rechargeable storage device. The Lit – ion device uses the movement of lithium ions between two electrodes, the positive and the negative. The Lit – ion batteries are very useful for consumer electronics and renewable grid storage. The Lit – ion batteries are known as high energy density device. They have cathode and anode ports. The anode port is constructed from graphite which stores the lithium ions and lets them to be free during discharge process. The cathode is usually made of lithium metal oxides (NMC or LFP) and functions as stores ions medium. Additional other two Lit – ion components are electrolyte and separator [1]. The electrolyte is a liquid or gel which carries the positively charged Lit – ions between the electrodes. The separator is a permeable physical barrier which keeps the electrodes apart. It prevents the short situation. By implementing charging process, the external electrical source forces the lithium ions to move back from the cathode to the anode. By implementing discharging process, the lithium ions move from the anode to the cathode. The movement is done through electrolyte by pushing electrons through an external circuit to power device. The advantages by using Lit – ion batteries are the high energy density, very low self – discharge rates, and no “memory effect”. If we compare Lit – ion batteries to nickel cadmium (NiCd) batteries, there is possibility to partially charge Lit – ion devices at any time without losing usable capacity or altering their voltage. The reason of Lit – ion batteries “not remember” is based on the fact that there is no crystal formation and better partial charging ability (keeping Lit – ion battery to 20% - 80% capacity extends its lifespan). The Lit – ion disadvantages are high cost, lifespan of two to three years, temperature sensitivity which cause to poor functionality under extreme cold or overheating (extreme heat). The Lit – ion batteries additionally have high risk of fire (thermal runaway). It happened when they mechanically damaged or over charged. Extreme heat or overheating can cause a short circuit which is call “thermal runaway”. It is uncontrolled chemical reaction which creates intense fires or explosions. The Lit – ion batteries are more expensive compare to NiCd batteries and have limited lifespan. They are very sensitive to temperature changes which effect their functionality. The important in controlling Lit – ion batteries to get conclusion on state of charge (SoC) and state of health (SoH) statuses. The lithium – ion battery state of charge (SoC) is the level of available charge relative to its maximum capacity [2]. It is ranging from zero (empty state) to one hundred percent (full state). It is battery fuel gauge and related to battery chemistry characteristics. We measure Lit – ion battery SoC by using battery management system (BMS) which run dedicated algorithms. The SoC measurement methods are coulomb counting (integrating the battery flowing current in and out), open circuit voltage OCV (measure the resting voltage of the battery since the Lit – ion batteries have flat voltage curve for discharge cycle), and combination of two methods alongside with temperature sensors to get the highest precision and the desire safety. The Lit – ion battery state of health (SoH) measure the device current electrochemical performance compare to its original current. We mark battery SoH by percentage which gives information on battery’s maximum capacity decay and increasing internal resistance for battery’s lifetime. Three regions characterized Lithium – ion batteries, the first is positive electrode, the second is separator, and the third is the negative electrode. The model simulates Lithium – ion battery in real time. A spaced discretized equations describe the real time of Lit – ion batteries. Optimization process is needed for the model. If we have a battery uniform current distribution then Lit – ion electrolyte concentration is expressed by partial differential equations [3]. The main variable is the electrolyte concentration and other important parameters are pore wall flux for negative and positive battery electrodes ( $j_n$ and $j_p$ , respectively). The battery pore wall flux is the rate which lithium ions move between the liquid electrolyte and the active solid material particles per unit area. It is boundary condition which related to coupling solid – phase diffusion and liquid phase electrolyte conservation equations. The source or sink of ions is determined by the pure well flux. It is implemented in both solid phases and liquid phases. The electrolyte concentration is the main Lit – ion variable when we want to inspect the battery operation dynamic. The electrolyte concentration variables are $C_{\text{ave}}^{\text{cathode}}$ , $C_{\text{ave}}$ , $C_{\text{ave}}^{\text{anode}}$ and presented by ordinary differential equations (ODEs) in time. The flux (concentration derivatives) of Lit – ion battery two interfaces are located at $x=l_p$ and at $x=l_p+l_s$ , respectively. The variable $x$ is related to Lit – ion sandwich distance from the first location at $x=0$ of positive electrode. We consider the battery concentration profile as a parabolic and the flux at the battery structure interfaces can be approximated. The Lit – ion battery solid phase diffusion approximations are reduced the computational efforts we need for real time BMS functionality. We get transient response of the model. The battery model helps us to predict the dynamical behaviour of the battery device in three regions which related to electrolyte concentration. The Lit – ion batteries discharge behaviour is the first outcome from battery model simulation. The Battery solid – phase diffusion coefficient ( $D_s$ ) dives us indication how fast lithium ions are moving inside electrode particles. If we define in the process delay parameters ( ${\tau }_i\forall i=1,2,3$ for ${\tau }_i\in ℝ\text{ ; }{\tau }_i>0$ ), they depend heavily on particle geometry, pulse duration, and relaxation period. Additionally, the Lit – ion electrode geometry modelling can define the internal time delay parameters which are diffusion time constants. It is for both microscopic and macroscopic scales. The defined model time delay parameters represent the actual physical time which required for Lit – ion battery lithium concentrations and electrical potentials to reorganize and get to steady – state (ss) equilibrium state during charge, discharge, or to the remain process periods. The Lit – ion battery electrolyte concentrations are delayed in time due to the solid – phase diffusion approximations and the transient response. The delayed in time parameters ( ${\tau }_i\forall i=1,2,3$ for ${\tau }_i\in ℝ\text{ ; }{\tau }_i>0$ ) Influence the Lit – ion battery system stability and need to be inspected deeply by stability criteria and non – linear dynamic theory. The standard local stability and stability switching are influenced by system overall parameters and model delay parameters.

Lithium-ion batteries solid-phase diffusion Approximations, transient response, electrolyte concentrations delayed in Time

We can define the Li – ion batteries electrolyte concentration as volume averaged over specific region, where C is the electrolyte concentration, we get the specific region concentration as $C_{\text{ave}}^{\text{cathode}}=\frac{1}{L_p}{\displaystyle \underset{x=0}{\overset{L_p}{\int }}Cdx}$ , $C_{\text{ave}}=\frac{1}{L_s}{\displaystyle \underset{L_p}{\overset{L_p+L_s}{\int }}C(x)dx}$  and $C_{\text{ave}}^{\text{anode}}=\frac{1}{L_n}{\displaystyle \underset{L_n+L_s}{\overset{L_p+L_s+L_n}{\int }}C(x)dx}$ , where $L_p=\frac{l_p}{L}$ , $L_s=\frac{l_s}{L}$ , $L_n=\frac{l_n}{L}$  , and $L=l_p+l_s+l_n$ . We have two interfaces in the Lit – ion cell sandwich structure, first is at $x=l_p$ and the second at $x=l_p+l_s$ . $L$ is the overall length of the Lit – ion positive and negative electrodes plus the separator length. $l_p$ is the length of positive electrode (thickness of the cathode), $l_s$ is the length of separator (thickness of the separator), and $l_n$ is the length of negative electrode (thickness of the anode). Typical values for $l_p$ , $l_s$  and $l_n$ are $80\times {\text{ 10}}^{-6}\text{ m}$ , $25\times {\text{ 10}}^{-6}\text{ m}$  and $88\times {\text{ 10}}^{-6}\text{ m}$ , respectively [4]. The concentration derivatives (flux) are at two interfaces, $x=l_p$ and at $x=l_p+l_s$ , respectively. We consider the Lit – ion concentration profile is parabolic and approximate the related fluxes. We can express the Lit – ion specific region concentrations by rate equations for getting the transient response of the Lit – ion model (Fig. 1). The related delay variables in time are $C_{\text{ave}}(t)\to C_{\text{ave}}(t-{\tau }_1)$ , $C_{\text{ave}}^{\text{cathode}}(t)\to C_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)$  and $C_{\text{ave}}^{\text{anode}}(t)\to C_{\text{ave}}^{\text{anode}}(t-{\tau }_3)$ , ${\tau }_1,{\tau }_2,{\tau }_3>0$ and ${\tau }_1,{\tau }_2,{\tau }_3\in ℝ$ .

 The Lit – ion rate delay differential equations (DDEs) are $\frac{d{C}_{\text{ave}}}{dt}=\frac{-6{\epsilon }_s^3(2C_{\text{ave}}(t-{\tau }_1)-C_1-C_3)}{L_s^2}$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3(C_1-C_{\text{ave}}^{\text{cathode}}(t-{\tau }_2))}{L_p^2}+J_p$ which related to positive electrode DDE for $C_{\text{ave}}^{\text{cathode}}$  , and $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3(C_3-C_{\text{ave}}^{\text{anode}}(t-{\tau }_3))(1-L_p-L_s)}{L_n^3}+J_n$  which related to negative electrode DDE for $C_{\text{ave}}^{\text{anode}}$ . We consider that the derivatives $\frac{d{C}_{\text{ave}}}{dt}$ , $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}$  and $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}$ do not affected by the delay parameters ( ${\tau }_1,{\tau }_2,{\tau }_3$ ). $J_n$ And $J_p$ are defined as the wall flux at the negative and positive electrodes, respectively [5]. In physics-based lithium – ion battery model ${\epsilon }_n$ and ${\epsilon }_p$ parameters are the volume fractions or porosities of the negative electrode (anode), and positive electrode (cathode), respectively. The negative Lit – ion battery ${\epsilon }_n$ parameter typical values are between 0.3 and 0.5 (30% to 50% porosity) values for negative electrode ( ${\epsilon }_n$ ). It is dependent on calendaring and whether it is a high – power or high energy cell. The positive Lit – ion battery ${\epsilon }_p$  parameter typical values are between 0.25 to 0.4 (25% to 40% porosity). We define $X$ parameter as $X=\frac{x}{L}$ , $J_p$ parameter as $J_p=[L^2\frac{(1-t_{+})}{c_{0}D{\epsilon }_p}]\times (\frac{I}{F{l}_p})$ , $J_n$ parameter as $J_n=[L^2\frac{(1-t_{+})}{c_{0}D{\epsilon }_n}]\times (\frac{I}{F{l}_n})$ . In Lithium – ion battery model, D parameter is the diffusion coefficient which is measured in $\frac{{\text{m}}^{\text{2}}}{\text{sec}}$ and represents how fast lithium ions move through the electrode particles or the electrolyte. It dictates the battery’s charge/discharge rate and power capacity. Typically, the D parameter value is $7.5\times {10}^{-10}$ (separator). In Lithium – ion battery model, F parameter is Faraday’s constant. It has typical value of $96487\frac{\text{C}}{\text{mol}}$ (separator) and represents the magnitude of electric charge per one mol of electrons. Lithium – ion F constant is related to ionic transport, diffusion and the model reaction kinetics. In lithium – ion battery model, C_O is the initial or reference concentration of lithium ions (measured by $\frac{\text{mol}}{{\text{m}}^{\text{3}}}$ ) [6]. The typical value of C_O  parameter is $1000\frac{\text{mol}}{{\text{m}}^{\text{3}}}$ (separator). It is also can be referred to battery solid phase or electrolyte phase. We define two functions $C_1$ and $C_3$  which are function of $C_{\text{ave}}^{\text{cathode}}$ , $C_{\text{ave}}$ $C_{\text{ave}}^{\text{anode}}$ variables and system overall parameters. I is defined as the applied battery cell current and it is related to the pore wall flux we get ( J or j ) [7]. The pore wall flux represents the local transfer rate of lithium ions between the solid active material and the liquid electrolyte. We define two functions ${\xi }_1={\xi }_1(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=C_1$ and ${\xi }_3={\xi }_3(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=C_3$ .

${\xi }_1={\xi }_1(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=C_1=\frac{3{\epsilon }_p^4{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)L_s+6{\epsilon }_s^4{L}_p{C}_{\text{ave}}(t-{\tau }_1)-2{\epsilon }_s^4{L}_p{C}_3}{3{\epsilon }_p^4{L}_s+4{\epsilon }_s^4{L}_p}\text{ (1)}$

${\xi }_1={\xi }_1(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=\frac{\begin{array}{l}3{\epsilon }_p^4{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)L_s+6{\epsilon }_s^4{L}_p{C}_{\text{ave}}(t-{\tau }_1)\\ -2{\epsilon }_s^4{L}_p{\xi }_3(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})\end{array}}{3{\epsilon }_p^4{L}_s+4{\epsilon }_s^4{L}_p}\text{ (2)}$

We define for simplicity some global system parameters: ${\Xi }_1=3{\epsilon }_p^4{L}_s$ ; ${\Xi }_2=6{\epsilon }_s^4{L}_p$ ; ${\Xi }_3=-2{\epsilon }_s^4{L}_p$ and

${\Xi }_4=3{\epsilon }_p^4{L}_s+4{\epsilon }_s^4{L}_p$ . Then we get the reduced ${\xi }_1(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})$ function:

${\xi }_1={\xi }_1(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=\frac{{\Xi }_1{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Xi }_2{C}_{\text{ave}}(t-{\tau }_1)+{\Xi }_3{\xi }_3(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})}{{\Xi }_4}\text{ (3)}$

${\xi }_3={\xi }_3(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=C_3=\frac{\begin{array}{l}-2{\epsilon }_s^4{L}_n{\epsilon }_p^4{L}_s{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+4{\epsilon }_s^8{L}_n{L}_p{C}_{\text{ave}}(t-{\tau }_1)+6{\epsilon }_s^4{L}_n{L}_s{\epsilon }_p^4{C}_{\text{ave}}(t-{\tau }_1)\\ +3{\epsilon }_n^4{L}_s^2{\epsilon }_p^4{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)+4{\epsilon }_s^4{L}_s{\epsilon }_s^4{L}_p{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)\end{array}}{4{\epsilon }_s^8{L}_n{L}_p+4{\epsilon }_s^4{L}_n{L}_s{\epsilon }_p^4+3{\epsilon }_n^4{L}_s^2{\epsilon }_p^4+4{\epsilon }_n^4{L}_s{\epsilon }_s^4{L}_p}\text{ (4)}$

${\xi }_3={\xi }_3(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=C_3=\frac{\begin{array}{l}-2{\epsilon }_s^4{L}_n{\epsilon }_p^4{L}_s{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+[4{\epsilon }_s^8{L}_n{L}_p+6{\epsilon }_s^4{L}_n{L}_s{\epsilon }_p^4]C_{\text{ave}}(t-{\tau }_1)\\ +[3{\epsilon }_n^4{L}_s^2{\epsilon }_p^4+4{\epsilon }_s^4{L}_s{\epsilon }_s^4{L}_p]C_{\text{ave}}^{\text{anode}}(t-{\tau }_3)\end{array}}{4{\epsilon }_s^8{L}_n{L}_p+4{\epsilon }_s^4{L}_n{L}_s{\epsilon }_p^4+3{\epsilon }_n^4{L}_s^2{\epsilon }_p^4+4{\epsilon }_n^4{L}_s{\epsilon }_s^4{L}_p}\text{ (5)}$

We define for simplicity some global system parameters: ${\Xi }_5=-2{\epsilon }_s^4{L}_n{\epsilon }_p^4{L}_s$ ; ${\Xi }_6=4{\epsilon }_s^8{L}_n{L}_p+6{\epsilon }_s^4{L}_n{L}_s{\epsilon }_p^4$ ,

${\Xi }_7=3{\epsilon }_n^4{L}_s^2{\epsilon }_p^4+4{\epsilon }_s^4{L}_s{\epsilon }_s^4{L}_p$ ; ${\Xi }_8=4{\epsilon }_s^8{L}_n{L}_p+4{\epsilon }_s^4{L}_n{L}_s{\epsilon }_p^4+3{\epsilon }_n^4{L}_s^2{\epsilon }_p^4+4{\epsilon }_n^4{L}_s{\epsilon }_s^4{L}_p$ .

${\xi }_3={\xi }_3(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})=C_3=\frac{{\Xi }_5{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Xi }_6{C}_{\text{ave}}(t-{\tau }_1)+{\Xi }_7{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)}{{\Xi }_8}\text{ (6)}$

We can summery our system DDEs and additional related functions:

$\frac{d{C}_{\text{ave}}}{dt}=\frac{-6{\epsilon }_s^3(2C_{\text{ave}}(t-{\tau }_1)-{\xi }_1-{\xi }_3)}{L_s^2}$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3({\xi }_1-C_{\text{ave}}^{\text{cathode}}(t-{\tau }_2))}{L_p^2}+J_p$ and

$\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3({\xi }_3-C_{\text{ave}}^{\text{anode}}(t-{\tau }_3))(1-L_p-L_s)}{L_n^3}+J_n$ ; ${\xi }_1=\frac{{\Xi }_1{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Xi }_2{C}_{\text{ave}}(t-{\tau }_1)+{\Xi }_3{\xi }_3}{{\Xi }_4}$

${\xi }_3=\frac{{\Xi }_5{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Xi }_6{C}_{\text{ave}}(t-{\tau }_1)+{\Xi }_7{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)}{{\Xi }_8}$

${\xi }_3=\frac{{\Xi }_5}{{\Xi }_8}{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+\frac{{\Xi }_6}{{\Xi }_8}{C}_{\text{ave}}(t-{\tau }_1)+\frac{{\Xi }_7}{{\Xi }_8}{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)$

We define for simplicity additional global parameter: ${\Upsilon }_1=\frac{{\Xi }_5}{{\Xi }_8}$ ; ${\Upsilon }_2=\frac{{\Xi }_6}{{\Xi }_8}$ ; ${\Upsilon }_3=\frac{{\Xi }_7}{{\Xi }_8}$ then

${\xi }_3={\Upsilon }_1{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Upsilon }_2{C}_{\text{ave}}(t-{\tau }_1)+{\Upsilon }_3{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)\text{ (7)}$

${\xi }_1=\frac{{\Xi }_1{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Xi }_2{C}_{\text{ave}}(t-{\tau }_1)+{\Xi }_3{\xi }_3}{{\Xi }_4}=\frac{\begin{array}{l}{\Xi }_1{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Xi }_2{C}_{\text{ave}}(t-{\tau }_1)\\ +{\Xi }_3[\frac{{\Xi }_5{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Xi }_6{C}_{\text{ave}}(t-{\tau }_1)+{\Xi }_7{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)}{{\Xi }_8}]\end{array}}{{\Xi }_4}\text{ (8)}$

${\xi }_1=\frac{[{\Xi }_1+\frac{{\Xi }_3{\Xi }_5}{{\Xi }_8}]C_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+[{\Xi }_2+\frac{{\Xi }_3{\Xi }_6}{{\Xi }_8}]C_{\text{ave}}(t-{\tau }_1)+\frac{{\Xi }_3{\Xi }_7}{{\Xi }_8}{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)}{{\Xi }_4}\text{ (9)}$

${\xi }_1=\frac{1}{{\Xi }_4}[{\Xi }_1+\frac{{\Xi }_3{\Xi }_5}{{\Xi }_8}]C_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+\frac{1}{{\Xi }_4}[{\Xi }_2+\frac{{\Xi }_3{\Xi }_6}{{\Xi }_8}]C_{\text{ave}}(t-{\tau }_1)+\frac{{\Xi }_3{\Xi }_7}{{\Xi }_8{\Xi }_4}{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)\text{ (10)}$

We define for simplicity additional global parameter: ${\Upsilon }_4=\frac{1}{{\Xi }_4}[{\Xi }_1+\frac{{\Xi }_3{\Xi }_5}{{\Xi }_8}]$ ; ${\Upsilon }_5=\frac{1}{{\Xi }_4}[{\Xi }_2+\frac{{\Xi }_3{\Xi }_6}{{\Xi }_8}]$ and

${\Upsilon }_6=\frac{{\Xi }_3{\Xi }_7}{{\Xi }_8{\Xi }_4}$ , then ${\xi }_1={\Upsilon }_4{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Upsilon }_5{C}_{\text{ave}}(t-{\tau }_1)+{\Upsilon }_6{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)$ .

We can summary our ${\xi }_1(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})$ and ${\xi }_3(C_{\text{ave}}^{\text{cathode}},C_{\text{ave}},C_{\text{ave}}^{\text{anode}})$ equations:

${\xi }_3={\Upsilon }_1{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Upsilon }_2{C}_{\text{ave}}(t-{\tau }_1)+{\Upsilon }_3{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)\text{ (11)}$

${\xi }_1={\Upsilon }_4{C}_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)+{\Upsilon }_5{C}_{\text{ave}}(t-{\tau }_1)+{\Upsilon }_6{C}_{\text{ave}}^{\text{anode}}(t-{\tau }_3)\text{ (12)}$

At fixed points (system equilibrium points) $\frac{d{C}_{\text{ave}}}{dt}=0$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=0$ ; $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=0$ , $C_{\text{ave}}(t-{\tau }_1)\approx C_{\text{ave}}(t)$

$\forall t\gg {\tau }_1\text{ ; }{\tau }_1\in ℝ$ , $C_{\text{ave}}^{\text{cathode}}(t-{\tau }_2)\approx C_{\text{ave}}^{\text{cathode}}(t)\forall t\gg {\tau }_2\text{ ; }{\tau }_2\in ℝ$ , and $C_{\text{ave}}^{\text{anode}}(t-{\tau }_3)\approx C_{\text{ave}}^{\text{anode}}(t)$

$\forall t\gg {\tau }_3\text{ ; }{\tau }_3\in ℝ$ ; ${\tau }_1,{\tau }_2,{\tau }_3>0$ . At fixed points ${\xi }_3^{*}={\Upsilon }_1[C_{\text{ave}}^{\text{cathode}}]*+{\Upsilon }_2[C_{\text{ave}}]*+{\Upsilon }_3[C_{\text{ave}}^{\text{anode}}]*$ and

${\xi }_1^{*}={\Upsilon }_4{[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_5{[C_{\text{ave}}]}^{*}+{\Upsilon }_6{[C_{\text{ave}}^{\text{anode}}]}^{*}$

$\frac{d{C}_{\text{ave}}}{dt}=0\Rightarrow \frac{-6{\epsilon }_s^3(2[C_{\text{ave}}]*-{\xi }_1^{*}-{\xi }_3^{*})}{L_s^2}=0$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=0\Rightarrow \frac{3{\epsilon }_p^3({\xi }_1^{*}-{[C_{\text{ave}}^{\text{cathode}}]}^{*})}{L_p^2}+J_p=0$ (13)

$\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=0\Rightarrow \frac{3{\epsilon }_n^3({\xi }_3^{*}-{[C_{\text{ave}}^{\text{anode}}]}^{*})(1-L_p-L_s)}{L_n^3}+J_n=0\text{ (14)}$

First fixed point equation (equilibrium points):

$\frac{-6{\epsilon }_s^3(2[C_{\text{ave}}]*-{\xi }_1^{*}-{\xi }_3^{*})}{L_s^2}=0\Rightarrow \frac{\begin{array}{l}-6{\epsilon }_s^3(2[C_{\text{ave}}]*-\{{\Upsilon }_4{[}^{C_{\text{ave}}^{\text{cathode}}}+{\Upsilon }_5{[}^{C_{\text{ave}}}+{\Upsilon }_6{[}^{C_{\text{ave}}^{\text{anode}}}\}\\ -\{{\Upsilon }_1[C_{\text{ave}}^{\text{cathode}}]*+{\Upsilon }_2[C_{\text{ave}}]*+{\Upsilon }_3[C_{\text{ave}}^{\text{anode}}]*\})\end{array}}{L_s^2}=0\text{ (15)}$

$\frac{-6{\epsilon }_s^3(2[C_{\text{ave}}]*-{\xi }_1^{*}-{\xi }_3^{*})}{L_s^2}=0\Rightarrow \frac{\begin{array}{l}-6{\epsilon }_s^3\{(2-{\Upsilon }_5-{\Upsilon }_2)[C_{\text{ave}}]*-({\Upsilon }_4+{\Upsilon }_1)[C_{\text{ave}}^{\text{cathode}}]*\\ -({\Upsilon }_6+{\Upsilon }_3)[C_{\text{ave}}^{\text{anode}}]*\}\end{array}}{L_s^2}=0\text{ (16)}$

$L_s^2\ne 0\text{ ; }-6{\epsilon }_s^3\ne 0\text{ ; }(2-{\Upsilon }_5-{\Upsilon }_2)[C_{\text{ave}}]*-({\Upsilon }_4+{\Upsilon }_1)[C_{\text{ave}}^{\text{cathode}}]*-({\Upsilon }_6+{\Upsilon }_3)[C_{\text{ave}}^{\text{anode}}]*=0\text{ (17)}$

Second fixed point equation (equilibrium points):

$\frac{3{\epsilon }_p^3({\xi }_1^{*}-{[C_{\text{ave}}^{\text{cathode}}]}^{*})}{L_p^2}+J_p=0\Rightarrow \frac{\begin{array}{l}3{\epsilon }_p^3(\{{\Upsilon }_4{[}^{C_{\text{ave}}^{\text{cathode}}}+{\Upsilon }_5{[}^{C_{\text{ave}}}+{\Upsilon }_6{[}^{C_{\text{ave}}^{\text{anode}}}\}\\ -{[}^{C_{\text{ave}}^{\text{cathode}}})\end{array}}{L_p^2}+J_p=0\text{ (18)}$

$\frac{3{\epsilon }_p^3({\xi }_1^{*}-{[C_{\text{ave}}^{\text{cathode}}]}^{*})}{L_p^2}+J_p=0\Rightarrow \frac{3{\epsilon }_p^3\{({\Upsilon }_4-1){[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_5{[C_{\text{ave}}]}^{*}+{\Upsilon }_6{[C_{\text{ave}}^{\text{anode}}]}^{*}\}}{L_p^2}+J_p=0\text{ (19)}$

$({\Upsilon }_4-1){[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_5{[C_{\text{ave}}]}^{*}+{\Upsilon }_6{[C_{\text{ave}}^{\text{anode}}]}^{*}+\frac{J_p{L}_p^2}{3{\epsilon }_p^3}=0\text{ (20)}$

Third fixed point equation (equilibrium points):

$\frac{3{\epsilon }_n^3({\xi }_3^{*}-{[C_{\text{ave}}^{\text{anode}}]}^{*})(1-L_p-L_s)}{L_n^3}+J_n=0\Rightarrow \frac{\begin{array}{l}3{\epsilon }_n^3(\{{\Upsilon }_1[C_{\text{ave}}^{\text{cathode}}]*+{\Upsilon }_2[C_{\text{ave}}]*\\ +{\Upsilon }_3[C_{\text{ave}}^{\text{anode}}]*\}-{[}^{C_{\text{ave}}^{\text{anode}}})(1-L_p-L_s)\end{array}}{L_n^3}+J_n=0\text{ (21)}$

$\frac{3{\epsilon }_n^3({\xi }_3^{*}-{[C_{\text{ave}}^{\text{anode}}]}^{*})(1-L_p-L_s)}{L_n^3}+J_n=0\Rightarrow \frac{\begin{array}{l}3{\epsilon }_n^3\{{\Upsilon }_1[C_{\text{ave}}^{\text{cathode}}]*+{\Upsilon }_2[C_{\text{ave}}]*\\ +({\Upsilon }_3-1){[}^{C_{\text{ave}}^{\text{anode}}}\}(1-L_p-L_s)\end{array}}{L_n^3}+J_n=0\text{ (22)}$

${\Upsilon }_1[C_{\text{ave}}^{\text{cathode}}]*+{\Upsilon }_2[C_{\text{ave}}]*+({\Upsilon }_3-1)[C_{\text{ave}}^{\text{anode}}]*+\frac{J_n{L}_n^3}{3{\epsilon }_n^3(1-L_p-L_s)}=0\text{ (23)}$

We can summery our system three fixed points equations:

$(2-{\Upsilon }_5-{\Upsilon }_2)[C_{\text{ave}}]*-({\Upsilon }_4+{\Upsilon }_1)[C_{\text{ave}}^{\text{cathode}}]*-({\Upsilon }_6+{\Upsilon }_3)[C_{\text{ave}}^{\text{anode}}]*=0\text{ (24)}$

$({\Upsilon }_4-1){[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_5{[C_{\text{ave}}]}^{*}+{\Upsilon }_6{[C_{\text{ave}}^{\text{anode}}]}^{*}+\frac{J_p{L}_p^2}{3{\epsilon }_p^3}=0\text{ (25)}$

${\Upsilon }_1[C_{\text{ave}}^{\text{cathode}}]*+{\Upsilon }_2[C_{\text{ave}}]*+({\Upsilon }_3-1)[C_{\text{ave}}^{\text{anode}}]*+\frac{J_n{L}_n^3}{3{\epsilon }_n^3(1-L_p-L_s)}=0\text{ (26)}$

We get system fixed points: $E^{*}({[C_{\text{ave}}]}^{*},{[C_{\text{ave}}^{\text{cathode}}]}^{*},[C_{\text{ave}}^{\text{anode}}]*)$ .

The Lit – ion battery delay differential equations (DDEs) model can be compared to the traditional reduced order Lit – ion battery model. The Lit – ion reduced – order model (ROM) simplifies the model’s partial differential equations (PDEs) in terms of solid diffusion and electrolyte migration. It is done into lower – dimensional approximations approach. It is applicable to single particle model (SPM) and SPM with electrode. We can implement it also to Galerkin projection methods. The Galerkin projection method is implemented as a mathematical technique in our case and use to convert continuous differential equations into a smaller and define – dimensional system equations. It is done by projecting the residual orthogonally elements onto a subspace of test function. Compare to Lit – ion reduced order model (ROM) the delay differential equations (DDEs) Lit – ion battery model uses time – delayed states such as diffusion lags, currents, or past voltages to account for thermal incubation, transport time or delayed internal feedback loops.

Lithium-ion batteries solid-phase diffusion Approximations, transient response, electrolyte concentrations characteristic equation τ1=τ, τ2=τ3=0

We get the expressions for ${\xi }_1={\xi }_1(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)$ and ${\xi }_3={\xi }_3(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)$ , then

${\xi }_3(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\Upsilon }_1{C}_{\text{ave}}^{\text{cathode}}(t)+{\Upsilon }_2{C}_{\text{ave}}(t-\tau )+{\Upsilon }_3{C}_{\text{ave}}^{\text{anode}}(t)\text{ (27)}$

${\xi }_1(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\Upsilon }_4{C}_{\text{ave}}^{\text{cathode}}(t)+{\Upsilon }_5{C}_{\text{ave}}(t-\tau )+{\Upsilon }_6{C}_{\text{ave}}^{\text{anode}}(t)\text{ (28)}$

The standard local stability analysis about any one of the equilibrium points of Lithium-ion batteries solid-phase electrolyte concentrations system consists in adding to coordinate [ $C_{\text{ave}}$ , $C_{\text{ave}}^{\text{cathode}}$ , $C_{\text{ave}}^{\text{anode}}$ ] arbitrarily small increments of exponential form $[C_{\text{ave}}{C}_{\text{ave}}^{\text{cathode}}{C}_{\text{ave}}^{\text{anode}}]e^{\lambda t}$ and retaining the first order term in $C_{\text{ave}}{C}_{\text{ave}}^{\text{cathode}}{C}_{\text{ave}}^{\text{anode}}$ . The system of three homogeneous equations leads to a polynomial characteristic equation in the eigenvalue $\lambda$ [8][9]. The Lithium-ion batteries solid-phase electrolyte concentrations system fixed values with arbitrarily small increments of exponential form $[C_{\text{ave}}{C}_{\text{ave}}^{\text{cathode}}{C}_{\text{ave}}^{\text{anode}}]e^{\lambda t}$ are $k=1$ (first fixed point), $k=2$ (second fixed point), $k=3$ (third fixed point), ...

$C_{\text{ave}}(t)={[C_{\text{ave}}]}^{*}+c_{\text{ave}}{e}^{\lambda t}$ ; $C_{\text{ave}}^{\text{cathode}}(t)=[C_{\text{ave}}^{\text{cathode}}]*+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}$ ; $C_{\text{ave}}^{\text{anode}}(t)=[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda t}$ (29)

$C_{\text{ave}}(t-\tau )={[C_{\text{ave}}]}^{*}+c_{\text{ave}}{e}^{\lambda (t-\tau )}$ , $C_{\text{ave}}^{\text{cathode}}(t-\tau )={[C_{\text{ave}}^{\text{cathode}}]}^{*}+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda (t-\tau )}$ , $C_{\text{ave}}^{\text{anode}}(t-\tau )=[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}$ (30)

$\begin{array}{c}{\xi }_3\left(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0\right)={\text{Υ}}_1\left({\left[C_{ave }^{cathode }\right]}^{\text{*}}+c_{ave }^{cathode }{e}^{\lambda t}\right)+{\text{Υ}}_2\left({\left[C_{ave }\right]}^{\text{*}}+c_{ave }{e}^{\lambda (t-\tau )}\right)\\ +{\text{Υ}}_3\left({\left[C_{ave }^{anode }\right]}^{\text{*}}+c_{ave }^{anode }{e}^{\lambda t}\right)\end{array}(31)$

$\begin{array}{l}{\xi }_3(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\Upsilon }_1{[}^{C_{\text{ave}}^{\text{cathode}}}+{\Upsilon }_2{[}^{C_{\text{ave}}}+{\Upsilon }_3[C_{\text{ave}}^{\text{anode}}]*\text{ (32)}\\ +{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}\end{array}$

We define ${\xi }_3^{*}$ fixed expression as ${\xi }_3^{*}={\Upsilon }_1{[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_2{[C_{\text{ave}}]}^{*}+{\Upsilon }_3[C_{\text{ave}}^{\text{anode}}]*$ then ${\xi }_3\left(t,{\tau }_1=\tau ,{\tau }_2={\tau }_3=0\right)$ is ${\xi }_3(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}$ .

$\begin{array}{l}{\xi }_1(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\Upsilon }_4({[}^{C_{\text{ave}}^{\text{cathode}}}+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t})+{\Upsilon }_5({[}^{C_{\text{ave}}}+c_{\text{ave}}{e}^{\lambda (t-\tau )})\text{ (33)}\\ +{\Upsilon }_6([C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda t})\end{array}$

$\begin{array}{l}{\xi }_1(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\Upsilon }_4{[}^{C_{\text{ave}}^{\text{cathode}}}+{\Upsilon }_5{[}^{C_{\text{ave}}}+{\Upsilon }_6[C_{\text{ave}}^{\text{anode}}]*\text{ (34)}\\ +{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}\end{array}$

We define ${\xi }_1^{*}$ fixed expression as ${\xi }_1^{*}={\Upsilon }_4{[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_5{[C_{\text{ave}}]}^{*}+{\Upsilon }_6[C_{\text{ave}}^{\text{anode}}]*$ then ${\xi }_1(t,{\tau }_1=\tau ,{\tau }_2={\tau }_3=0)$

is ${\xi }_1(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}$ .

We choose the above expressions for our $C_{\text{ave}}^{\text{cathode}}(t),C_{\text{ave}}(t),C_{\text{ave}}^{\text{anode}}(t)$ as small displacement $c_{\text{ave}},c_{\text{ave}}^{\text{cathode}},c_{\text{ave}}^{\text{anode}}$ from the system fixed points at time $t=0$ ( $C_{\text{ave}}=C_{\text{ave}}(t)$ , $C_{\text{ave}}^{\text{cathode}}=C_{\text{ave}}^{\text{cathode}}(t)$ , $C_{\text{ave}}^{\text{anode}}=C_{\text{ave}}^{\text{anode}}(t)$ ).

$\frac{d{C}_{\text{ave}}(t)}{dt}=c_{\text{ave}}\lambda e^{\lambda t}$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}(t)}{dt}=c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}$ ; $\frac{d{C}_{\text{ave}}^{\text{anode}}(t)}{dt}=c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}$ (35)

$\frac{d{C}_{\text{ave}}(t-\tau )}{dt}=c_{\text{ave}}\lambda e^{\lambda (t-\tau )}$ , $\frac{d{C}_{\text{ave}}^{\text{cathode}}(t-\tau )}{dt}=c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda (t-\tau )}$ , $\frac{d{C}_{\text{ave}}^{\text{anode}}(t-\tau )}{dt}=c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda (t-\tau )}$ (36)

Under our specific delay parameters value ( ${\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0$ ), we get the following system delay differential equations (DDEs):

$\frac{d{C}_{\text{ave}}}{dt}=\frac{-6{\epsilon }_s^3(2C_{\text{ave}}(t-\tau )-{\xi }_1-{\xi }_3)}{L_s^2}$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3({\xi }_1-C_{\text{ave}}^{\text{cathode}}(t))}{L_p^2}+J_p$ and $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3({\xi }_3-C_{\text{ave}}^{\text{anode}}(t))(1-L_p-L_s)}{L_n^3}+J_n$ . The above ${\xi }_1$ , ${\xi }_3$ expressions are related to our specific case ( $\left({\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0\right),{\xi }_1\left(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0\right)$ ), and ${\xi }_3(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)$ .

$\frac{d{C}_{\text{ave}}}{dt}=\frac{-6{\epsilon }_s^3(2C_{\text{ave}}(t-\tau )-{\xi }_1-{\xi }_3)}{L_s^2}\Rightarrow c_{\text{ave}}\lambda e^{\lambda t}=\frac{\begin{array}{l}-6{\epsilon }_s^3\{2({[}^{C_{\text{ave}}}+c_{\text{ave}}{e}^{\lambda (t-\tau )})\\ -({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t})\\ -({\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t})\}\end{array}}{L_s^2}\text{ (37)}$

$\begin{array}{l}-6{\epsilon }_s^3\{2({[}^{C_{\text{ave}}}+c_{\text{ave}}{e}^{\lambda (t-\tau )})-({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t})\\ -({\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t})\}=-6{\epsilon }_s^{3}2{[}^{C_{\text{ave}}}-6{\epsilon }_s^{3}2c_{\text{ave}}{e}^{\lambda (t-\tau )}\text{ (38)}\\ +6{\epsilon }_s^3{\xi }_1^{*}+6{\epsilon }_s^3{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+6{\epsilon }_s^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}\\ +6{\epsilon }_s^3{\xi }_3^{*}+6{\epsilon }_s^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+6{\epsilon }_s^3{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}\end{array}$

$\begin{array}{l}-6{\epsilon }_s^{3}2{[}^{C_{\text{ave}}}-6{\epsilon }_s^{3}2c_{\text{ave}}{e}^{\lambda (t-\tau )}+6{\epsilon }_s^3{\xi }_1^{*}+6{\epsilon }_s^3{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}\\ +6{\epsilon }_s^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}+6{\epsilon }_s^3{\xi }_3^{*}+6{\epsilon }_s^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+6{\epsilon }_s^3{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}\text{ (39)}\\ =-6({\epsilon }_s^{3}2{[}^{C_{\text{ave}}}-{\epsilon }_s^3{\xi }_1^{*}-{\epsilon }_s^3{\xi }_3^{*})+6\{{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]c_{\text{ave}}{e}^{-\lambda \tau }\\ +{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]c_{\text{ave}}^{\text{cathode}}+{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}\}{e}^{\lambda t}\end{array}$

We get the following first eigenvalue function:

$c_{\text{ave}}\lambda e^{\lambda t}=\frac{\begin{array}{l}-6({\epsilon }_s^{3}2{[}^{C_{\text{ave}}}-{\epsilon }_s^3{\xi }_1^{*}-{\epsilon }_s^3{\xi }_3^{*})+6\{{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]c_{\text{ave}}{e}^{-\lambda \tau }\\ +{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]c_{\text{ave}}^{\text{cathode}}+{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}\}{e}^{\lambda t}\end{array}}{L_s^2}\text{ (40)}$

$c_{\text{ave}}\lambda e^{\lambda t}=\frac{-6{\epsilon }_s^3(2{[C_{\text{ave}}]}^{*}-{\xi }_1^{*}-{\xi }_3^{*})}{L_s^2}+\frac{\begin{array}{l}6\{{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]c_{\text{ave}}{e}^{-\lambda \tau }\\ +{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]c_{\text{ave}}^{\text{cathode}}+{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}\}{e}^{\lambda t}\end{array}}{L_s^2}\text{ (41)}$

Since $\frac{d{C}_{ave }}{dt}=\frac{-6{\epsilon }_s^3\left(2C_{ave }(t-\tau )-{\xi }_1-{\xi }_3\right)}{L_s^2}$ , at fixed point ( $\frac{d{C}_{\text{ave}}}{dt}=0$ ; $\underset{\begin{array}{l}t\to \infty \\ t\gg \tau \end{array}}{\lim }{C}_{\text{ave}}(t-\tau )=C_{\text{ave}}(t)$ ) exists $\frac{-6{\epsilon }_s^3(2{[C_{\text{ave}}]}^{*}-{\xi }_1^{*}-{\xi }_3^{*})}{L_s^2}=0$ and the first eigenvalue equation is

$c_{\text{ave}}\lambda =\frac{6\{{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]c_{\text{ave}}{e}^{-\lambda \tau }+{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]c_{\text{ave}}^{\text{cathode}}+{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}\}}{L_s^2}\text{ (42)}$

$(\frac{6{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]e^{-\lambda \tau }}{L_s^2}-\lambda )c_{\text{ave}}+\frac{6{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]}{L_s^2}{c}_{\text{ave}}^{\text{cathode}}+\frac{6{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]}{L_s^2}{c}_{\text{ave}}^{\text{anode}}=0\text{ (43)}$

We define for simplicity system global parameters: ${\Xi }_1=\frac{6{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]}{L_s^2}$ ; ${\Xi }_2=\frac{6{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]}{L_s^2}$

${\Xi }_3=\frac{6{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]}{L_s^2}$ then the first eigenvalue equation is $({\Xi }_1{e}^{-\lambda \tau }-\lambda )c_{\text{ave}}+{\Xi }_2{c}_{\text{ave}}^{\text{cathode}}+{\Xi }_3{c}_{\text{ave}}^{\text{anode}}=0$ .

$\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3({\xi }_1-C_{\text{ave}}^{\text{cathode}}(t))}{L_p^2}+J_p\Rightarrow c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_p^3[({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}\\ +{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t})-({[}^{C_{\text{ave}}^{\text{cathode}}}+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t})]\end{array}}{L_p^2}+J_p\text{ (44)}$

$\begin{array}{l}3{\epsilon }_p^3[({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t})-({[}^{C_{\text{ave}}^{\text{cathode}}}+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t})]=3{\epsilon }_p^3{\xi }_1^{*}\text{ (45)}\\ +3{\epsilon }_p^3{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}-3{\epsilon }_p^3{[}^{C_{\text{ave}}^{\text{cathode}}}-3{\epsilon }_p^3{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}\end{array}$

$\begin{array}{l}3{\epsilon }_p^3{\xi }_1^{*}+3{\epsilon }_p^3{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}-3{\epsilon }_p^3{[}^{C_{\text{ave}}^{\text{cathode}}}-3{\epsilon }_p^3{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}\text{ (46)}\\ =3{\epsilon }_p^3({\xi }_1^{*}-{[}^{C_{\text{ave}}^{\text{cathode}}})+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}\end{array}$

We get the following second eigenvalue function:

$c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_p^3({\xi }_1^{*}-{[}^{C_{\text{ave}}^{\text{cathode}}})+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}\\ +3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}\end{array}}{L_p^2}+J_p\text{ (47)}$

$c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}=\frac{3{\epsilon }_p^3({\xi }_1^{*}-{[C_{\text{ave}}^{\text{cathode}}]}^{*})}{L_p^2}+J_p+\frac{3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}}{L_p^2}\text{ (48)}$

Since $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3({\xi }_1-C_{\text{ave}}^{\text{cathode}}(t))}{L_p^2}+J_p$ , at fixed point ( $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=0$ ) exists

$\frac{3{\epsilon }_p^3({\xi }_1^{*}-{[C_{\text{ave}}^{\text{cathode}}]}^{*})}{L_p^2}+J_p=0$ and the second eigenvalue equation is

$c_{\text{ave}}^{\text{cathode}}\lambda =\frac{3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{-\lambda \tau }+3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}}{L_p^2}\text{ (49)}$

$\frac{3{\epsilon }_p^3{\Upsilon }_5{e}^{-\lambda \tau }}{L_p^2}{c}_{\text{ave}}+(\frac{3{\epsilon }_p^3({\Upsilon }_4-1)}{L_p^2}-\lambda )c_{\text{ave}}^{\text{cathode}}+\frac{3{\epsilon }_p^3{\Upsilon }_6}{L_p^2}{c}_{\text{ave}}^{\text{anode}}=0\text{ (50)}$

We define for simplicity system global parameters: ${\Xi }_4=\frac{3{\epsilon }_p^3{\Upsilon }_5}{L_p^2}$ ; ${\Xi }_5=\frac{3{\epsilon }_p^3({\Upsilon }_4-1)}{L_p^2}$ ; ${\Xi }_6=\frac{3{\epsilon }_p^3{\Upsilon }_6}{L_p^2}$ then the second eigenvalue equation is ${\text{Ξ}}_4{e}^{-\lambda \tau }{c}_{ave }+\left({\text{Ξ}}_5-\lambda \right)c_{ave }^{cathode }+{\text{Ξ}}_6{c}_{ave }^{anode }=0$ .

$\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3({\xi }_3-C_{\text{ave}}^{\text{anode}}(t))(1-L_p-L_s)}{L_n^3}+J_n\Rightarrow c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_n^3([{\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}]\\ -[[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda t}])(1-L_p-L_s)\end{array}}{L_n^3}+J_n\text{ (51)}$

$\begin{array}{l}3{\epsilon }_n^3([{\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}]-[[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda t}])(1-L_p-L_s)\hfill \\ \begin{array}{l}=3{\epsilon }_n^3{\xi }_3^{*}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\\ +3{\epsilon }_n^3{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}(1-L_p-L_s)-3{\epsilon }_n^3[C_{\text{ave}}^{\text{anode}}]*(1-L_p-L_s)-3{\epsilon }_n^3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}(1-L_p-L_s)\end{array}\hfill \end{array}\text{ (52)}$

$\begin{array}{l}3{\epsilon }_n^3{\xi }_3^{*}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\hfill \\ \begin{array}{l}+3{\epsilon }_n^3{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}(1-L_p-L_s)-3{\epsilon }_n^3[C_{\text{ave}}^{\text{anode}}]*(1-L_p-L_s)-3{\epsilon }_n^3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda t}(1-L_p-L_s)\\ =3{\epsilon }_n^3({\xi }_3^{*}-[C_{\text{ave}}^{\text{anode}}]*)(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\\ +3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3[{\Upsilon }_3-1]c_{\text{ave}}^{\text{anode}}{e}^{\lambda t}(1-L_p-L_s)\end{array}\hfill \end{array}\text{ (53)}$

We get the following third eigenvalue function:

$c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_n^3({\xi }_3^{*}-[C_{\text{ave}}^{\text{anode}}]*)(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\\ +3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3[{\Upsilon }_3-1]c_{\text{ave}}^{\text{anode}}{e}^{\lambda t}(1-L_p-L_s)\end{array}}{L_n^3}+J_n\text{ (54)}$

$\begin{array}{l}{c}_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}=\frac{3{\epsilon }_n^3({\xi }_3^{*}-[C_{\text{ave}}^{\text{anode}}]*)(1-L_p-L_s)}{L_n^3}+J_n+\frac{3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)}{L_n^3}\text{ (55)}\\ +\frac{3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)}{L_n^3}+\frac{3{\epsilon }_n^3[{\Upsilon }_3-1]c_{\text{ave}}^{\text{anode}}{e}^{\lambda t}(1-L_p-L_s)}{L_n^3}\end{array}$

Since $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3({\xi }_3-C_{\text{ave}}^{\text{anode}}(t))(1-L_p-L_s)}{L_n^3}+J_n$ , at fixed point ( $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=0$ ) exists

$\frac{3{\epsilon }_n^3\left({\xi }_3^{\text{*}}-{\left[C_{ave }^{anode }\right]}^{\text{*}}\right)\left(1-L_p-L_s\right)}{L_n^3}+J_n=0$ and the third eigenvalue equation [10][11] is

$c_{\text{ave}}^{\text{anode}}\lambda =\frac{3{\epsilon }_n^3{\Upsilon }_2(1-L_p-L_s)}{L_n^3}{e}^{-\lambda \tau }{c}_{\text{ave}}+\frac{3{\epsilon }_n^3{\Upsilon }_1(1-L_p-L_s)}{L_n^3}{c}_{\text{ave}}^{\text{cathode}}+\frac{3{\epsilon }_n^3[{\Upsilon }_3-1](1-L_p-L_s)}{L_n^3}{c}_{\text{ave}}^{\text{anode}}\text{ (56)}$

We define for simplicity system global parameters: ${\Xi }_7=\frac{3{\epsilon }_n^3{\Upsilon }_2(1-L_p-L_s)}{L_n^3}$ ; ${\Xi }_8=\frac{3{\epsilon }_n^3{\Upsilon }_1(1-L_p-L_s)}{L_n^3}$

${\Xi }_9=\frac{3{\epsilon }_n^3[{\Upsilon }_3-1](1-L_p-L_s)}{L_n^3}$ then the third eigenvalue equation is ${\text{Ξ}}_7{e}^{-\lambda \tau }{c}_{ave }+{\text{Ξ}}_8{c}_{ave }^{cathode }+\left({\text{Ξ}}_9-\lambda \right)c_{ave }^{anode }=0$ .

We can summary our system three eigenvalues equation with delay parameter ( $\tau$ ):

First eigenvalue equation: $({\Xi }_1{e}^{-\lambda \tau }-\lambda )c_{\text{ave}}+{\Xi }_2{c}_{\text{ave}}^{\text{cathode}}+{\Xi }_3{c}_{\text{ave}}^{\text{anode}}=0$

Second eigenvalue equation: ${\text{Ξ}}_4{e}^{-\lambda \tau }{c}_{ave }+\left({\text{Ξ}}_5-\lambda \right)c_{ave }^{cathode }+{\text{Ξ}}_6{c}_{ave }^{anode }=0$

Third eigenvalue equation: ${\text{Ξ}}_7{e}^{-\lambda \tau }{c}_{ave }+{\text{Ξ}}_8{c}_{ave }^{cathode }+\left({\text{Ξ}}_9-\lambda \right)c_{ave }^{anode }=0$

The eigenvalue matrix we get is $\left(\begin{array}{ccc}\left({\text{Ξ}}_1{e}^{-\lambda \tau }-\lambda \right)& {\text{Ξ}}_2& {\text{Ξ}}_3\\ {\text{Ξ}}_4{e}^{-\lambda \tau }& \left({\text{Ξ}}_5-\lambda \right)& {\text{Ξ}}_6\\ {\text{Ξ}}_7{e}^{-\lambda \tau }& {\text{Ξ}}_8& \left({\text{Ξ}}_9-\lambda \right)\end{array}\right)\left(\begin{array}{l}{c}_{ave }\hfill \\ c_{ave }^{cathode }\hfill \\ c_{ave }^{anode }\hfill \end{array}\right)=0$

$A-\lambda I=\left(\begin{array}{ccc}({\Xi }_1{e}^{-\lambda \tau }-\lambda )& {\Xi }_2& {\Xi }_3\\ {\Xi }_4{e}^{-\lambda \tau }& ({\Xi }_5-\lambda )& {\Xi }_6\\ {\Xi }_7{e}^{-\lambda \tau }& {\Xi }_8& ({\Xi }_9-\lambda )\end{array}\right)$ ; $\det (A-\lambda I)=0$ (57)

$\begin{array}{l}\det (A-\lambda I)=({\Xi }_1{e}^{-\lambda \tau }-\lambda )\det \left(\begin{array}{cc}({\Xi }_5-\lambda )& {\Xi }_6\\ {\Xi }_8& ({\Xi }_9-\lambda )\end{array}\right)-{\Xi }_2\det \left(\begin{array}{cc}{\Xi }_4{e}^{-\lambda \tau }& {\Xi }_6\\ {\Xi }_7{e}^{-\lambda \tau }& ({\Xi }_9-\lambda )\end{array}\right)\text{ (58)}\\ +{\Xi }_3\det \left(\begin{array}{cc}{\Xi }_4{e}^{-\lambda \tau }& ({\Xi }_5-\lambda )\\ {\Xi }_7{e}^{-\lambda \tau }& {\Xi }_8\end{array}\right)\end{array}$

$\begin{array}{l}\det (A-\lambda I)=({\Xi }_1{e}^{-\lambda \tau }-\lambda )[({\Xi }_5-\lambda )({\Xi }_9-\lambda )-{\Xi }_8{\Xi }_6]-{\Xi }_2[{\Xi }_4({\Xi }_9-\lambda )e^{-\lambda \tau }-{\Xi }_7{\Xi }_6{e}^{-\lambda \tau }]\text{ (59)}\\ +{\Xi }_3[{\Xi }_8{\Xi }_4{e}^{-\lambda \tau }-{\Xi }_7({\Xi }_5-\lambda )e^{-\lambda \tau }]\end{array}$

$\begin{array}{l}\det (A-\lambda I)=({\Xi }_1{e}^{-\lambda \tau }-\lambda )([{\Xi }_5{\Xi }_9-{\Xi }_8{\Xi }_6]-({\Xi }_5+{\Xi }_9)\lambda +{\lambda }^2)-{\Xi }_2({\Xi }_4{\Xi }_9{e}^{-\lambda \tau }-{\Xi }_4\lambda e^{-\lambda \tau }-{\Xi }_7{\Xi }_6{e}^{-\lambda \tau })\text{ (60)}\\ +{\Xi }_3({\Xi }_8{\Xi }_4{e}^{-\lambda \tau }-{\Xi }_7{\Xi }_5{e}^{-\lambda \tau }+{\Xi }_7\lambda e^{-\lambda \tau })\end{array}$

$\begin{array}{l}\det (A-\lambda I)=[{\Xi }_5{\Xi }_9-{\Xi }_8{\Xi }_6]{\Xi }_1{e}^{-\lambda \tau }-({\Xi }_5+{\Xi }_9)\lambda {\Xi }_1{e}^{-\lambda \tau }+{\lambda }^2{\Xi }_1{e}^{-\lambda \tau }-[{\Xi }_5{\Xi }_9-{\Xi }_8{\Xi }_6]\lambda \\ +({\Xi }_5+{\Xi }_9){\lambda }^2-{\lambda }^3-{\Xi }_2{\Xi }_4{\Xi }_9{e}^{-\lambda \tau }+{\Xi }_2{\Xi }_4\lambda e^{-\lambda \tau }+{\Xi }_2{\Xi }_7{\Xi }_6{e}^{-\lambda \tau }+{\Xi }_3{\Xi }_8{\Xi }_4{e}^{-\lambda \tau }-{\Xi }_3{\Xi }_7{\Xi }_5{e}^{-\lambda \tau }\text{ (61)}\\ +{\Xi }_3{\Xi }_7\lambda e^{-\lambda \tau }\end{array}$

$\begin{array}{l}\det (A-\lambda I)=-[{\Xi }_5{\Xi }_9-{\Xi }_8{\Xi }_6]\lambda +({\Xi }_5+{\Xi }_9){\lambda }^2-{\lambda }^3\\ +\{[{\Xi }_5{\Xi }_9-{\Xi }_8{\Xi }_6]{\Xi }_1-{\Xi }_2{\Xi }_4{\Xi }_9+{\Xi }_2{\Xi }_7{\Xi }_6+{\Xi }_3{\Xi }_8{\Xi }_4-{\Xi }_3{\Xi }_7{\Xi }_5\text{ (62)}\\ +({\Xi }_2{\Xi }_4-({\Xi }_5+{\Xi }_9){\Xi }_1+{\Xi }_3{\Xi }_7)\lambda +{\lambda }^2{\Xi }_1\}{e}^{-\lambda \tau }\end{array}$

We define for simplicity global parameters: ${\Omega }_1=-[{\Xi }_5{\Xi }_9-{\Xi }_8{\Xi }_6]$ ; ${\Omega }_2={\Xi }_5+{\Xi }_9$ ; ${\Omega }_4={\Xi }_2{\Xi }_4-({\Xi }_5+{\Xi }_9){\Xi }_1+{\Xi }_3{\Xi }_7$ ,

${\Omega }_3=[{\Xi }_5{\Xi }_9-{\Xi }_8{\Xi }_6]{\Xi }_1-{\Xi }_2{\Xi }_4{\Xi }_9+{\Xi }_2{\Xi }_7{\Xi }_6+{\Xi }_3{\Xi }_8{\Xi }_4-{\Xi }_3{\Xi }_7{\Xi }_5$

$\det (A-\lambda I)={\Omega }_1\lambda +{\Omega }_2{\lambda }^2-{\lambda }^3+[{\Omega }_3+{\Omega }_4\lambda +{\Xi }_1{\lambda }^2]e^{-\lambda \tau }$ ; $\det (A-\lambda I)=0$ (63)

The general characteristics equation is $D(\lambda ,\tau )={\Omega }_1\lambda +{\Omega }_2{\lambda }^2-{\lambda }^3+[{\Omega }_3+{\Omega }_4\lambda +{\Xi }_1{\lambda }^2]e^{-\lambda \tau }$ We study the occurrence of any possible stability switching resulting from the increase of value of the time delay for the general characteristics equation $D(\lambda ,\tau )$ [12][13]. The characteristics equation is $D(\lambda ,\tau )=P_n(\lambda ,\tau )+Q_m(\lambda ,\tau )e^{-\lambda \tau }$ ; the expression for $P_n(\lambda ,\tau )$ : $P_n(\lambda ,\tau )={\displaystyle \sum }_{k=0}^n{\lambda }^k{p}_k(\tau )$ . The expression for $Q_m(\lambda ,\tau )$ : $Q_m(\lambda ,\tau )={\displaystyle \sum }_{k=0}^m{\lambda }^k{q}_k(\tau )$ . The characteristics equation is $D(\lambda ,\tau )=P_n(\lambda ,\tau )+Q_m(\lambda ,\tau )e^{-\lambda \tau }$ where $n=3$ , $m=2$ and $n>m$ . The expression for $P_n(\lambda ,\tau )$ : $P_n(\lambda ,\tau )={\displaystyle \sum }_{k=0}^{n=3}{\lambda }^k{p}_k(\tau )$ ; $p_0(\tau )=0$ ; $p_1(\tau )={\Omega }_1$ , $p_2(\tau )={\Omega }_2$ ; $p_3(\tau )=-1$ . The expression for $Q_m(\lambda ,\tau )$ : $Q_m(\lambda ,\tau )={\displaystyle \sum }_{k=0}^{m=2}{\lambda }^k{q}_k(\tau )$ ; $q_0(\tau )={\Omega }_3$ ; $q_1(\tau )={\Omega }_4$ , $q_2(\tau )={\Xi }_1$ . The homogeneous system for $C_{\text{ave}}$ , $C_{\text{ave}}^{\text{cathode}}$ , $C_{\text{ave}}^{\text{anode}}$ leads to a characteristic equation for the eigenvalue $\lambda$ . (1) If $\lambda =j\omega$ , $\omega \in ℝ$ then $P(j\omega )+Q(j\omega )\ne 0$ , (2) $\left|\frac{Q(\lambda )}{P(\lambda )}\right|$ is bounded for $\left|\lambda \right|\to \infty$ , $\mathrm{Re}\lambda \ge 0$ and no roots bifurcation from $\infty$ , (3) $F(\omega )=\left|P(j\omega ){|}^2-\right|Q(j\omega ){|}^2$ has a finite number of zeros. Indeed, this is a polynomial in $\omega$ , (4) Each $\omega (q_i,q_k)$ of $F(\omega )=0$ is continuous and differentiable with respect to $q_i$ , $q_k$ .

$P_n(\lambda =j\omega ,\tau )={\displaystyle \sum }_{k=0}^{n=3}{\lambda }^k{p}_k(\tau )=-{\Omega }_2{\omega }^2+j({\omega }^3+{\Omega }_1\omega )$ ; $Q_m(\lambda =j\omega ,\tau )={\displaystyle \sum }_{k=0}^{m=2}{\lambda }^k{q}_k(\tau )={\Omega }_3-{\Xi }_1{\omega }^2+j{\Omega }_4\omega$ (64)

$\begin{array}{c}�F(\omega )=|P(j\omega ){|}^2-|Q(j\omega ){|}^2={\left(-{\text{Ω}}_2{\omega }^2\right)}^2+{\left({\omega }^3+{\text{Ω}}_1\omega \right)}^2-\left[{\left({\text{Ω}}_3-{\text{Ξ}}_1{\omega }^2\right)}^2+{\left({\text{Ω}}_4\omega \right)}^2\right]={\omega }^6(65)\\ � +\left(2{\text{Ω}}_1+{\text{Ω}}_2^2\right){\omega }^4+{\text{Ω}}_1^2{\omega }^2-\left[{\text{Ω}}_3^2-2{\text{Ω}}_3{\text{Ξ}}_1{\omega }^2+{\text{Ξ}}_1^2{\omega }^4+{\text{Ω}}_4^2{\omega }^2\right]\end{array}$

$F(\omega )=\left|P(j\omega ){|}^2-\right|Q(j\omega ){|}^2={\omega }^6+(2{\Omega }_1+{\Omega }_2^2){\omega }^4+{\Omega }_1^2{\omega }^2-{\Omega }_3^2+2{\Omega }_3{\Xi }_1{\omega }^2-{\Xi }_1^2{\omega }^4-{\Omega }_4^2{\omega }^2\text{ (66)}$

$F(\omega )=\left|P(j\omega ){|}^2-\right|Q(j\omega ){|}^2={\omega }^6+[(2{\Omega }_1+{\Omega }_2^2)-{\Xi }_1^2]{\omega }^4+[{\Omega }_1^2+2{\Omega }_3{\Xi }_1-{\Omega }_4^2]{\omega }^2-{\Omega }_3^2\text{ (67)}$

We define ${\Phi }_0=-{\Omega }_3^2$ ; ${\Phi }_2={\Omega }_1^2+2{\Omega }_3{\Xi }_1-{\Omega }_4^2$ ; ${\Phi }_4=(2{\Omega }_1+{\Omega }_2^2)-{\Xi }_1^2$ ; ${\Phi }_6=1$ . $F(\omega )=0$ implies ${\displaystyle \sum }_{k=0}^3{\Phi }_{2k}{\omega }^{2k}=0$ . Its roots are given by solving the polynomial ( ${\displaystyle \sum }_{k=0}^3{\Phi }_{2k}{\omega }^{2k}=0$ ). $P_R(j\omega ,\tau )=-{\Omega }_2{\omega }^2$ , $P_I(j\omega ,\tau )={\omega }^3+{\Omega }_1\omega$ , $Q_R(j\omega ,\tau )={\Omega }_3-{\Xi }_1{\omega }^2$ , $Q_I(j\omega ,\tau )={\Omega }_4\omega$ .

$\sin \theta (\tau )=\sin \theta (\omega \tau )=g(\omega )$ ; $\cos \theta (\tau )=\cos \theta (\omega \tau )=h(\omega )$ (68)

$\sin \theta (\tau )=\frac{-P_R(j\omega ,\tau )Q_I(j\omega ,\tau )+P_I(j\omega ,\tau )Q_R(j\omega ,\tau )}{|Q(j\omega ,\tau ){|}^2}\text{ (69)}$

$\cos \theta (\tau )=-\frac{[P_R(j\omega ,\tau )Q_R(j\omega ,\tau )+P_I(j\omega ,\tau )Q_I(j\omega ,\tau )]}{|Q(j\omega ,\tau ){|}^2}\text{ (70)}$

The above $\sin \theta (\tau )$ and $\cos \theta (\tau )$ expressions are jointly with $F(\omega )=0\Rightarrow {\displaystyle \sum }_{k=0}^3{\Phi }_{2k}{\omega }^{2k}=0$ . Eigenvalue $\lambda =0$ is not a root of the characteristic equation. Furthermore, $P(\lambda )$ , $Q(\lambda )$ are analytic functions as the coefficients in $P$ and $Q$ are real. We have $\overline{P(-j\omega )}=P(j\omega )$ and $\overline{Q(-j\omega )}=Q(j\omega )$ thus $\lambda =j\omega ;\omega >0$ may be an eigenvalue of the characteristic equation [14][15]. The analysis consists of identifying the roots of the characteristic equation situated on the imaginary axis of the complex $\lambda$ plane, by increasing the parameters ${\epsilon }_s,L_s,{\epsilon }_p,L_p,{\epsilon }_n,L_n,J_p,J_n,\dots ,\text{Re}\lambda$ $\mathrm{Re}\lambda$ may, at the crossing, change its sign from (-) to (+), that is from a stable focus E* to an unstable, or vice versa. This feature may be further assessed by examining the sign of the partial derivatives with respect to ${\epsilon }_s,L_s,{\epsilon }_p,L_p,{\epsilon }_n,L_n,J_p,J_n$ parameters. Upon separating into real and imaginary parts, with $P=P_R+j{P}_I$ ; $Q=Q_R+j{Q}_I$ and $P_{\omega }=P_{R\omega }+j{P}_{I\omega }$ ; $Q_{\omega }=Q_{R\omega }+j{Q}_{I\omega }$ ; $P^2=P_R^2+P_I^2$ .

When ( $x$ ) can be any Lithium-ion batteries solid-phase electrolyte concentrations system parameters ${\epsilon }_s,L_s,{\epsilon }_p,L_p,{\epsilon }_n,L_n,J_p,J_n,\dots$ , ... and time delay parameter $\tau$ .

$F_{\omega }=2[(P_{R\omega }{P}_R+P_{I\omega }{P}_I)-(Q_{R\omega }{Q}_R+Q_{I\omega }{Q}_I)]$ ; $F_x=2[(P_{Rx}{P}_R+P_{Ix}{P}_I)-(Q_{Rx}{Q}_R+Q_{Ix}{Q}_I)]$ ; $\begin{array}{c}{\omega }_x=-\frac{F_x}{F_{\omega }}\end{array}$ (71)

We choose our specific parameter as time delay $x=\tau$ then $P_R=-{\Omega }_2{\omega }^2$ ; $P_I={\omega }^3+{\Omega }_1\omega$ ; $Q_R={\Omega }_3-{\Xi }_1{\omega }^2$ , $Q_I={\Omega }_4\omega$ . The related derivatives are $P_{R\omega }=-2{\Omega }_2\omega$ ; $P_{I\omega }=3{\omega }^2+{\Omega }_1$ ; $Q_{R\omega }=-2{\Xi }_1\omega$ , $Q_{I\omega }={\Omega }_4$ .

$F_{\tau }=2[(P_{R\tau }{P}_R+P_{I\tau }{P}_I)-(Q_{R\tau }{Q}_R+Q_{I\tau }{Q}_I)]=0$ ; $V=\left(P_R{P}_{I\tau }-P_I{P}_{R\tau }\right)-\left(Q_R{Q}_{I\tau }-Q_I{Q}_{R\tau }\right)=0$ (72)

$U=(P_R{P}_{I\omega }-P_I{P}_{R\omega })-(Q_R{Q}_{I\omega }-Q_I{Q}_{R\omega })\text{ (73)}$

$U=(-{\Omega }_2{\omega }^2[3{\omega }^2+{\Omega }_1]+2{\Omega }_2{\omega }^2[{\omega }^2+{\Omega }_1])-({\Omega }_4[{\Omega }_3-{\Xi }_1{\omega }^2]+2{\Omega }_4{\Xi }_1{\omega }^2)\text{ (74)}$

$U=-{\Omega }_2{\omega }^2[3{\omega }^2+{\Omega }_1]+2{\Omega }_2{\omega }^2[{\omega }^2+{\Omega }_1]-{\Omega }_4[{\Omega }_3-{\Xi }_1{\omega }^2]-2{\Omega }_4{\Xi }_1{\omega }^2\text{ (75)}$

By differentiating $F(\omega ,\tau )$ with respect to $\tau$ , we get $F_{\omega }\frac{\partial \omega }{\partial \tau }+F_{\tau }=0$ ; $\tau \in I\Rightarrow \frac{\partial \omega }{\partial \tau }=-\frac{F_{\tau }}{F_{\omega }}$ .

${\Lambda }^{-1}(\tau )={(\frac{\partial \mathrm{Re}\lambda }{\partial \tau })}_{\lambda =j\omega }$ ; $\frac{\partial \omega }{\partial \tau }={\omega }_{\tau }=-\frac{F_{\tau }}{F_{\omega }}$ ; ${\Lambda }^{-1}(\tau )=\mathrm{Re}\left\{\frac{-2[U+\tau |P{|}^2]+j{F}_{\omega }}{F_{\tau }+j2[V+\omega |P{|}^2]}\right\}$ (76)

$\text{sign}\{{\Lambda }^{-1}(\tau )\}=\text{sign}\left\{F_{\omega }\right\}\text{sign}\{\tau \frac{\partial \omega }{\partial \tau }+\omega +\frac{U\frac{\partial \omega }{\partial \tau }+V}{|P{|}^2}\}\text{ (77)}$

Since $F_{\tau }=0\Rightarrow \frac{\partial \omega }{\partial \tau }={\omega }_{\tau }=-\frac{F_{\tau }}{F_{\omega }}{|}_{F_{\omega }\ne 0}=0$ and $V=0$ then $\text{sign}\{{\Lambda }^{-1}(\tau )\}=\text{sign}\left\{F_{\omega }\right\}\text{sign}\left\{\omega \right\}$ . The following ( $\text{sign}\left\{{\Lambda }^{-1}(\tau )\right\}$ ) examines the possibility of stability transitions (bifurcations) Lithium-ion batteries solid-phase electrolyte concentrations system, about the equilibrium point E* as a result of a variation of delay $\tau$ . We identify the roots of our system characteristic equation situated on the imaginary axis of the complex

$\lambda$ - plane whereby increasing the delay parameter $\tau$ , $\text{Re}\lambda$ may at the crossing, change its sign from – to +, that is from a stable focus E* to an unstable one, or vice versa. We get the expression for $F(\omega ,\tau )$ Lithium-ion batteries solid-phase electrolyte concentrations system parameter values. We define those $\omega ,\tau$ values which fulfill $F(\omega ,\tau )=0$ . We ignore negative, complex, and imaginary values of $\omega$ for specific $\tau$ values [16][17].

Lithium-ion batteries solid-phase diffusion Approximations, transient response, electrolyte concentrations characteristic equation τ1=τ2=0, τ3=τ

We get the expressions for ${\xi }_1={\xi }_1(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )$ and ${\xi }_3={\xi }_3(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )$ , then

${\xi }_3={\Upsilon }_1{C}_{\text{ave}}^{\text{cathode}}(t)+{\Upsilon }_2{C}_{\text{ave}}(t)+{\Upsilon }_3{C}_{\text{ave}}^{\text{anode}}(t-\tau )\text{ (78)}$

${\xi }_1={\Upsilon }_4{C}_{\text{ave}}^{\text{cathode}}(t)+{\Upsilon }_5{C}_{\text{ave}}(t)+{\Upsilon }_6{C}_{\text{ave}}^{\text{anode}}(t-\tau )\text{ (79)}$

The standard local stability analysis about any one of the equilibrium points of Lithium-ion batteries solid-phase electrolyte concentrations system consists in adding to coordinate [ $C_{\text{ave}}$ , $C_{\text{ave}}^{\text{cathode}}$ , $C_{\text{ave}}^{\text{anode}}$ ] arbitrarily small increments of exponential form $\left[c_{ave },c_{ave }^{cathode },c_{ave }^{anode }\right]e^{\lambda t}$ and retaining the first order term in $c_{ave },c_{ave }^{cathode },c_{ave }^{anode }$ . The system of three homogeneous equations leads to a polynomial characteristic equation in the eigenvalue $\lambda$ [18][19]. The Lithium-ion batteries solid-phase electrolyte concentrations system fixed values with arbitrarily small increments of exponential form $\left[c_{ave },c_{ave }^{cathode },c_{ave }^{anode }\right]e^{\lambda t}$ are $k=1$ (first fixed point), $k=2$ (second fixed point), $k=3$ (third fixed point), ...

$C_{\text{ave}}(t)={[C_{\text{ave}}]}^{*}+c_{\text{ave}}{e}^{\lambda t}$ ; $C_{\text{ave}}^{\text{cathode}}(t)=[C_{\text{ave}}^{\text{cathode}}]*+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}$ ; $C_{\text{ave}}^{\text{anode}}(t)=[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda t}$ (80)

$C_{\text{ave}}(t-\tau )={[C_{\text{ave}}]}^{*}+c_{\text{ave}}{e}^{\lambda (t-\tau )}$ , $C_{\text{ave}}^{\text{cathode}}(t-\tau )={[C_{\text{ave}}^{\text{cathode}}]}^{*}+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda (t-\tau )}$ , $C_{\text{ave}}^{\text{anode}}(t-\tau )=[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}$ (81)

${\xi }_3={\Upsilon }_1([C_{\text{ave}}^{\text{cathode}}]*+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t})+{\Upsilon }_2({[C_{\text{ave}}]}^{*}+c_{\text{ave}}{e}^{\lambda t})+{\Upsilon }_3([C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\text{ (82)}$

$\begin{array}{l}{\xi }_3(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )={\Upsilon }_1{[}^{C_{\text{ave}}^{\text{cathode}}}+{\Upsilon }_2{[}^{C_{\text{ave}}}+{\Upsilon }_3[C_{\text{ave}}^{\text{anode}}]*\text{ (83)}\\ +{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}\end{array}$

We define ${\xi }_3^{*}$ fixed expression as ${\xi }_3^{*}={\Upsilon }_1{[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_2{[C_{\text{ave}}]}^{*}+{\Upsilon }_3[C_{\text{ave}}^{\text{anode}}]*$ then ${\xi }_3(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )$ is ${\xi }_3(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )={\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}$ .

$\begin{array}{l}{\xi }_1(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )={\Upsilon }_4([C_{\text{ave}}^{\text{cathode}}]*+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t})+{\Upsilon }_5({[C_{\text{ave}}]}^{*}+c_{\text{ave}}{e}^{\lambda t})\hfill \\ +{\Upsilon }_6([C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\hfill \end{array}\text{ (84)}$

$\begin{array}{l}{\xi }_1(t,{\tau }_1=\tau ,{\tau }_2=0,{\tau }_3=0)={\Upsilon }_4{[}^{C_{\text{ave}}^{\text{cathode}}}+{\Upsilon }_5{[}^{C_{\text{ave}}}+{\Upsilon }_6[C_{\text{ave}}^{\text{anode}}]*\text{ (85)}\\ +{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}\end{array}$

We define ${\xi }_1^{*}$ fixed expression as ${\xi }_1^{*}={\Upsilon }_4{[C_{\text{ave}}^{\text{cathode}}]}^{*}+{\Upsilon }_5{[C_{\text{ave}}]}^{*}+{\Upsilon }_6[C_{\text{ave}}^{\text{anode}}]*$ then ${\xi }_1(t,{\tau }_1=0,{\tau }_2={\tau }_3=\tau )$

is ${\xi }_1(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )={\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}$ .

We choose the above expressions for our $C_{\text{ave}}^{\text{cathode}}(t),C_{\text{ave}}(t),C_{\text{ave}}^{\text{anode}}(t)$ as small displacement $c_{\text{ave}},c_{\text{ave}}^{\text{cathode}},c_{\text{ave}}^{\text{anode}}$ from the system fixed points at time $t=0$ ( $C_{\text{ave}}=C_{\text{ave}}(t)$ , $C_{\text{ave}}^{\text{cathode}}=C_{\text{ave}}^{\text{cathode}}(t)$ , $C_{\text{ave}}^{\text{anode}}=C_{\text{ave}}^{\text{anode}}(t)$ ).

$\frac{d{C}_{\text{ave}}(t)}{dt}=c_{\text{ave}}\lambda e^{\lambda t}$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}(t)}{dt}=c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}$ ; $\frac{d{C}_{\text{ave}}^{\text{anode}}(t)}{dt}=c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}$ (86)

$\frac{d{C}_{\text{ave}}(t-\tau )}{dt}=c_{\text{ave}}\lambda e^{\lambda (t-\tau )}$ , $\frac{d{C}_{\text{ave}}^{\text{cathode}}(t-\tau )}{dt}=c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda (t-\tau )}$ , $\frac{d{C}_{\text{ave}}^{\text{anode}}(t-\tau )}{dt}=c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda (t-\tau )}$ (87)

Under our specific delay parameters value ( ${\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau$ ), we get the following system delay differential equations (DDEs):

$\frac{d{C}_{\text{ave}}}{dt}=\frac{-6{\epsilon }_s^3(2C_{\text{ave}}(t)-{\xi }_1-{\xi }_3)}{L_s^2}$ ; $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3({\xi }_1-C_{\text{ave}}^{\text{cathode}}(t))}{L_p^2}+J_p$ and $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3({\xi }_3-C_{\text{ave}}^{\text{anode}}(t-\tau ))(1-L_p-L_s)}{L_n^3}+J_n$ [20][21]. The above ${\xi }_1$ , ${\xi }_3$ expressions are related to our specific case ( ${\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau$ ), ${\xi }_1(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )$ and ${\xi }_3(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )$ .

$\frac{d{C}_{\text{ave}}}{dt}=\frac{-6{\epsilon }_s^3(2C_{\text{ave}}(t)-{\xi }_1-{\xi }_3)}{L_s^2}\Rightarrow c_{\text{ave}}\lambda e^{\lambda t}=\frac{\begin{array}{l}-6{\epsilon }_s^3\{2({[}^{C_{\text{ave}}}+c_{\text{ave}}{e}^{\lambda t})\\ -({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\\ -({\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\}\end{array}}{L_s^2}\text{ (88)}$

$\begin{array}{l}-6{\epsilon }_s^3\{2({[}^{C_{\text{ave}}}+c_{\text{ave}}{e}^{\lambda t})-({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\\ -({\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\}=-6{\epsilon }_s^{3}2{[}^{C_{\text{ave}}}-6{\epsilon }_s^{3}2c_{\text{ave}}{e}^{\lambda t}\text{ (89)}\\ +6{\epsilon }_s^3{\xi }_1^{*}+6{\epsilon }_s^3{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}\\ +6{\epsilon }_s^3{\xi }_3^{*}+6{\epsilon }_s^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+6{\epsilon }_s^3{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}\end{array}$

$\begin{array}{l}-6{\epsilon }_s^3\{2({[}^{C_{\text{ave}}}+c_{\text{ave}}{e}^{\lambda t})-({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\\ -({\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})\}=-6({\epsilon }_s^{3}2{[}^{C_{\text{ave}}}-{\epsilon }_s^3{\xi }_1^{*}-{\epsilon }_s^3{\xi }_3^{*})\text{ (90)}\\ +6\{{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]c_{\text{ave}}+{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]c_{\text{ave}}^{\text{cathode}}+{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }\}{e}^{\lambda t}\end{array}$

We get the following first eigenvalue function:

$c_{\text{ave}}\lambda e^{\lambda t}=\frac{\begin{array}{l}-6({\epsilon }_s^{3}2{[}^{C_{\text{ave}}}-{\epsilon }_s^3{\xi }_1^{*}-{\epsilon }_s^3{\xi }_3^{*})+6\{{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]c_{\text{ave}}\\ +{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]c_{\text{ave}}^{\text{cathode}}+{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }\}{e}^{\lambda t}\end{array}}{L_s^2}\text{ (91)}$

$c_{\text{ave}}\lambda e^{\lambda t}=\frac{-6{\epsilon }_s^3(2{[C_{\text{ave}}]}^{*}-{\xi }_1^{*}-{\xi }_3^{*})}{L_s^2}+\frac{\begin{array}{l}6\{{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]c_{\text{ave}}\\ +{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]c_{\text{ave}}^{\text{cathode}}+{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }\}{e}^{\lambda t}\end{array}}{L_s^2}\text{ (92)}$

Since $\frac{d{C}_{\text{ave}}}{dt}=\frac{-6{\epsilon }_s^3(2C_{\text{ave}}(t)-{\xi }_1-{\xi }_3)}{L_s^2}$ , at fixed point ( $\frac{d{C}_{\text{ave}}}{dt}=0$ ; $\underset{\begin{array}{l}t\to \infty \\ t\gg \tau \end{array}}{\lim }{C}_{\text{ave}}^{\text{anode}}(t-\tau )=C_{\text{ave}}^{\text{anode}}(t)$ ) exists $\frac{-6{\epsilon }_s^3(2{[C_{\text{ave}}]}^{*}-{\xi }_1^{*}-{\xi }_3^{*})}{L_s^2}=0$ and the first eigenvalue equation is

$c_{\text{ave}}\lambda =\frac{6{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]}{L_s^2}{c}_{\text{ave}}+\frac{6{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]}{L_s^2}{c}_{\text{ave}}^{\text{cathode}}+\frac{6{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }}{L_s^2}\text{ (93)}$

$(\frac{6{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]}{L_s^2}-\lambda )c_{\text{ave}}+\frac{6{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]}{L_s^2}{c}_{\text{ave}}^{\text{cathode}}+\frac{6{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]c_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }}{L_s^2}=0\text{ (94)}$

We define for simplicity system global parameters: ${\Xi }_1=\frac{6{\epsilon }_s^3[-2+{\Upsilon }_5+{\Upsilon }_2]}{L_s^2}$ ; ${\Xi }_2=\frac{6{\epsilon }_s^3[{\Upsilon }_4+{\Upsilon }_1]}{L_s^2}$

${\Xi }_3=\frac{6{\epsilon }_s^3[{\Upsilon }_6+{\Upsilon }_3]}{L_s^2}$ then the first eigenvalue equation is $({\Xi }_1-\lambda )c_{\text{ave}}+{\Xi }_2{c}_{\text{ave}}^{\text{cathode}}+{\Xi }_3{e}^{-\lambda \tau }{c}_{\text{ave}}^{\text{anode}}=0$ .

$\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3({\xi }_1-C_{\text{ave}}^{\text{cathode}}(t))}{L_p^2}+J_p\Rightarrow c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_p^3[({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}\\ +{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})-({[}^{C_{\text{ave}}^{\text{cathode}}}+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t})]\end{array}}{L_p^2}+J_p\text{ (95)}$

$\begin{array}{l}3{\epsilon }_p^3[({\xi }_1^{*}+{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )})-({[}^{C_{\text{ave}}^{\text{cathode}}}+c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t})]=3{\epsilon }_p^3{\xi }_1^{*}\text{ (96)}\\ +3{\epsilon }_p^3{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}-3{\epsilon }_p^3{[}^{C_{\text{ave}}^{\text{cathode}}}-3{\epsilon }_p^3{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}\end{array}$

$\begin{array}{l}3{\epsilon }_p^3{\xi }_1^{*}+3{\epsilon }_p^3{\Upsilon }_4{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}-3{\epsilon }_p^3{[}^{C_{\text{ave}}^{\text{cathode}}}-3{\epsilon }_p^3{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}\text{ (97)}\\ =3{\epsilon }_p^3({\xi }_1^{*}-{[}^{C_{\text{ave}}^{\text{cathode}}})+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}\end{array}$

We get the following second eigenvalue function:

$c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_p^3({\xi }_1^{*}-{[}^{C_{\text{ave}}^{\text{cathode}}})+3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}\\ +3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{c}_{\text{ave}}{e}^{\lambda (t-\tau )}\end{array}}{L_p^2}+J_p\text{ (98)}$

$c_{\text{ave}}^{\text{cathode}}\lambda e^{\lambda t}=\frac{3{\epsilon }_p^3({\xi }_1^{*}-{[C_{\text{ave}}^{\text{cathode}}]}^{*})}{L_p^2}+J_p+\frac{3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}{e}^{\lambda t}+3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}}{L_p^2}\text{ (99)}$

Since $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=\frac{3{\epsilon }_p^3({\xi }_1-C_{\text{ave}}^{\text{cathode}}(t))}{L_p^2}+J_p$ , at fixed point ( $\frac{d{C}_{\text{ave}}^{\text{cathode}}}{dt}=0$ ) exists

$\frac{3{\epsilon }_p^3({\xi }_1^{*}-{[C_{\text{ave}}^{\text{cathode}}]}^{*})}{L_p^2}+J_p=0$ and the second eigenvalue equation is

$c_{\text{ave}}^{\text{cathode}}\lambda =\frac{3{\epsilon }_p^3{\Upsilon }_5{c}_{\text{ave}}+3{\epsilon }_p^3({\Upsilon }_4-1)c_{\text{ave}}^{\text{cathode}}+3{\epsilon }_p^3{\Upsilon }_6{c}_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }}{L_p^2}\text{ (100)}$

$\frac{3{\epsilon }_p^3{\Upsilon }_5}{L_p^2}{c}_{\text{ave}}+(\frac{3{\epsilon }_p^3({\Upsilon }_4-1)}{L_p^2}-\lambda )c_{\text{ave}}^{\text{cathode}}+\frac{3{\epsilon }_p^3{\Upsilon }_6}{L_p^2}{c}_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }=0\text{ (101)}$

We define for simplicity system global parameters: ${\Xi }_4=\frac{3{\epsilon }_p^3{\Upsilon }_5}{L_p^2}$ ; ${\Xi }_5=\frac{3{\epsilon }_p^3({\Upsilon }_4-1)}{L_p^2}$ ; ${\Xi }_6=\frac{3{\epsilon }_p^3{\Upsilon }_6}{L_p^2}$ then the second eigenvalue equation is ${\Xi }_4{c}_{\text{ave}}+({\Xi }_5-\lambda )c_{\text{ave}}^{\text{cathode}}+{\Xi }_6{c}_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }=0$ .

$\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3({\xi }_3-C_{\text{ave}}^{\text{anode}}(t-\tau ))(1-L_p-L_s)}{L_n^3}+J_n\Rightarrow c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_n^3([{\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}]\\ -[[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}])(1-L_p-L_s)\end{array}}{L_n^3}+J_n\text{ (102)}$

Since ${\xi }_3(t,{\tau }_1=0,{\tau }_2=0,{\tau }_3=\tau )={\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}$

$\begin{array}{l}3{\epsilon }_n^3([{\xi }_3^{*}+{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}+{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}+{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}]-[[C_{\text{ave}}^{\text{anode}}]*+c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}])(1-L_p-L_s)\hfill \\ \begin{array}{l}=3{\epsilon }_n^3{\xi }_3^{*}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}(1-L_p-L_s)\\ +3{\epsilon }_n^3{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)-3{\epsilon }_n^3[C_{\text{ave}}^{\text{anode}}]*(1-L_p-L_s)-3{\epsilon }_n^3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\end{array}\hfill \end{array}\text{ (103)}$

$\begin{array}{l}3{\epsilon }_n^3{\xi }_3^{*}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}(1-L_p-L_s)\hfill \\ \begin{array}{l}+3{\epsilon }_n^3{\Upsilon }_3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)-3{\epsilon }_n^3[C_{\text{ave}}^{\text{anode}}]*(1-L_p-L_s)-3{\epsilon }_n^3{c}_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\\ =3{\epsilon }_n^3({\xi }_3^{*}-[C_{\text{ave}}^{\text{anode}}]*)(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}(1-L_p-L_s)\\ +3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3[{\Upsilon }_3-1]c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\end{array}\hfill \end{array}\text{ (104)}$

We get the following third eigenvalue function:

$c_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}=\frac{\begin{array}{l}3{\epsilon }_n^3({\xi }_3^{*}-[C_{\text{ave}}^{\text{anode}}]*)(1-L_p-L_s)+3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}(1-L_p-L_s)\\ +3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)+3{\epsilon }_n^3[{\Upsilon }_3-1]c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)\end{array}}{L_n^3}+J_n\text{ (105)}$

$\begin{array}{l}{c}_{\text{ave}}^{\text{anode}}\lambda e^{\lambda t}=\frac{3{\epsilon }_n^3({\xi }_3^{*}-[C_{\text{ave}}^{\text{anode}}]*)(1-L_p-L_s)}{L_n^3}+J_n+\frac{3{\epsilon }_n^3{\Upsilon }_2{c}_{\text{ave}}{e}^{\lambda t}(1-L_p-L_s)}{L_n^3}\text{ (106)}\\ +\frac{3{\epsilon }_n^3{\Upsilon }_1{c}_{\text{ave}}^{\text{cathode}}{e}^{\lambda t}(1-L_p-L_s)}{L_n^3}+\frac{3{\epsilon }_n^3[{\Upsilon }_3-1]c_{\text{ave}}^{\text{anode}}{e}^{\lambda (t-\tau )}(1-L_p-L_s)}{L_n^3}\end{array}$

Since $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=\frac{3{\epsilon }_n^3({\xi }_3-C_{\text{ave}}^{\text{anode}}(t))(1-L_p-L_s)}{L_n^3}+J_n$ , at fixed point ( $\frac{d{C}_{\text{ave}}^{\text{anode}}}{dt}=0$ ) exists

$\frac{3{\epsilon }_n^3({\xi }_3^{*}-[C_{\text{ave}}^{\text{anode}}]*)(1-L_p-L_s)}{L_n^3}=0$ and the third eigenvalue equation is

$c_{\text{ave}}^{\text{anode}}\lambda =\frac{3{\epsilon }_n^3{\Upsilon }_2(1-L_p-L_s)}{L_n^3}{c}_{\text{ave}}+\frac{3{\epsilon }_n^3{\Upsilon }_1(1-L_p-L_s)}{L_n^3}{c}_{\text{ave}}^{\text{cathode}}+\frac{3{\epsilon }_n^3[{\Upsilon }_3-1](1-L_p-L_s)}{L_n^3}{e}^{-\lambda \tau }{c}_{\text{ave}}^{\text{anode}}\text{ (107)}$

We define for simplicity system global parameters: ${\Xi }_7=\frac{3{\epsilon }_n^3{\Upsilon }_2(1-L_p-L_s)}{L_n^3}$ ; ${\Xi }_8=\frac{3{\epsilon }_n^3{\Upsilon }_1(1-L_p-L_s)}{L_n^3}$

${\Xi }_9=\frac{3{\epsilon }_n^3[{\Upsilon }_3-1](1-L_p-L_s)}{L_n^3}$ then the third eigenvalue equation is $c_{\text{ave}}^{\text{anode}}\lambda ={\Xi }_7{c}_{\text{ave}}+{\Xi }_8{c}_{\text{ave}}^{\text{cathode}}+{\Xi }_9{e}^{-\lambda \tau }{c}_{\text{ave}}^{\text{anode}}$ .

We can summary our system three eigenvalues equation [22][23] with delay parameter ( $\tau$ ):

First eigenvalue equation: $({\Xi }_1-\lambda )c_{\text{ave}}+{\Xi }_2{c}_{\text{ave}}^{\text{cathode}}+{\Xi }_3{e}^{-\lambda \tau }{c}_{\text{ave}}^{\text{anode}}=0$

Second eigenvalue equation: ${\Xi }_4{c}_{\text{ave}}+({\Xi }_5-\lambda )c_{\text{ave}}^{\text{cathode}}+{\Xi }_6{c}_{\text{ave}}^{\text{anode}}{e}^{-\lambda \tau }=0$

Third eigenvalue equation: $c_{\text{ave}}^{\text{anode}}\lambda ={\Xi }_7{c}_{\text{ave}}+{\Xi }_8{c}_{\text{ave}}^{\text{cathode}}+{\Xi }_9{e}^{-\lambda \tau }{c}_{\text{ave}}^{\text{anode}}$

The eigenvalue matrix we get is $\left(\begin{array}{ccc}({\Xi }_1-\lambda )& {\Xi }_2& {\Xi }_3{e}^{-\lambda \tau }\\ {\Xi }_4& ({\Xi }_5-\lambda )& {\Xi }_6{e}^{-\lambda \tau }\\ {\Xi }_7& {\Xi }_8& ({\Xi }_9-\lambda )e^{-\lambda \tau }\end{array}\right)\left(\begin{array}{l}{c}_{\text{ave}}\\ c_{\text{ave}}^{\text{cathode}}\\ c_{\text{ave}}^{\text{anode}}\end{array}\right)=0$

$A-\lambda I=\left(\begin{array}{ccc}({\Xi }_1-\lambda )& {\Xi }_2& {\Xi }_3{e}^{-\lambda \tau }\\ {\Xi }_4& ({\Xi }_5-\lambda )& {\Xi }_6{e}^{-\lambda \tau }\\ {\Xi }_7& {\Xi }_8& ({\Xi }_9-\lambda )e^{-\lambda \tau }\end{array}\right)$ ; $\det (A-\lambda I)=0$ (108)

$\begin{array}{l}\det (A-\lambda I)=({\Xi }_1-\lambda )\det \left(\begin{array}{cc}({\Xi }_5-\lambda )& {\Xi }_6{e}^{-\lambda \tau }\\ {\Xi }_8& ({\Xi }_9-\lambda )e^{-\lambda \tau }\end{array}\right)-{\Xi }_2\det \left(\begin{array}{cc}{\Xi }_4& {\Xi }_6{e}^{-\lambda \tau }\\ {\Xi }_7& ({\Xi }_9-\lambda )e^{-\lambda \tau }\end{array}\right)\text{ (109)}\\ +{\Xi }_3{e}^{-\lambda \tau }\det \left(\begin{array}{cc}{\Xi }_4& ({\Xi }_5-\lambda )\\ {\Xi }_7& {\Xi }_8\end{array}\right)\end{array}$

$\begin{array}{l}\det (A-\lambda I)=({\Xi }_1-\lambda )[({\Xi }_5-\lambda )({\Xi }_9-\lambda )e^{-\lambda \tau }-{\Xi }_8{\Xi }_6{e}^{-\lambda \tau }]-{\Xi }_2[{\Xi }_4({\Xi }_9-\lambda )e^{-\lambda \tau }-{\Xi }_7{\Xi }_6{e}^{-\lambda \tau }]\text{ (110)}\\ +{\Xi }_3{e}^{-\lambda \tau }[{\Xi }_4{\Xi }_8-{\Xi }_7({\Xi }_5-\lambda )]\end{array}$

$\begin{array}{l}\det (A-\lambda I)=({\Xi }_1-\lambda )[({\Xi }_5{\Xi }_9-\lambda ({\Xi }_5+{\Xi }_9)+{\lambda }^2)e^{-\lambda \tau }-{\Xi }_8{\Xi }_6{e}^{-\lambda \tau }]\\ -{\Xi }_2[{\Xi }_4{\Xi }_9{e}^{-\lambda \tau }-{\Xi }_4\lambda e^{-\lambda \tau }-{\Xi }_7{\Xi }_6{e}^{-\lambda \tau }]+{\Xi }_3{e}^{-\lambda \tau }[{\Xi }_4{\Xi }_8-{\Xi }_5{\Xi }_7+{\Xi }_7\lambda ]\text{ (111)}\end{array}$

$\begin{array}{l}\det (A-\lambda I)=({\Xi }_1{\Xi }_5{\Xi }_9-\lambda {\Xi }_1({\Xi }_5+{\Xi }_9)+{\Xi }_1{\lambda }^2)e^{-\lambda \tau }-{\Xi }_1{\Xi }_8{\Xi }_6{e}^{-\lambda \tau }\\ -\lambda ({\Xi }_5{\Xi }_9-\lambda ({\Xi }_5+{\Xi }_9)+{\lambda }^2)e^{-\lambda \tau }+{\Xi }_8{\Xi }_6\lambda e^{-\lambda \tau }-{\Xi }_2{\Xi }_4{\Xi }_9{e}^{-\lambda \tau }\\ +{\Xi }_2{\Xi }_4\lambda e^{-\lambda \tau }+{\Xi }_2{\Xi }_7{\Xi }_6{e}^{-\lambda \tau }+{\Xi }_4{\Xi }_8{\Xi }_3{e}^{-\lambda \tau }-{\Xi }_5{\Xi }_7{\Xi }_3{e}^{-\lambda \tau }+{\Xi }_7{\Xi }_3\lambda e^{-\lambda \tau }\text{ (112)}\end{array}$

$\begin{array}{l}\det (A-\lambda I)={\Xi }_1{\Xi }_5{\Xi }_9{e}^{-\lambda \tau }-{\Xi }_1({\Xi }_5+{\Xi }_9)\lambda e^{-\lambda \tau }+{\Xi }_1{\lambda }^2{e}^{-\lambda \tau }-{\Xi }_1{\Xi }_8{\Xi }_6{e}^{-\lambda \tau }\text{ (113)}\\ -{\Xi }_5{\Xi }_9\lambda e^{-\lambda \tau }+({\Xi }_5+{\Xi }_9){\lambda }^2{e}^{-\lambda \tau }-{\lambda }^3{e}^{-\lambda \tau }+{\Xi }_8{\Xi }_6\lambda e^{-\lambda \tau }-{\Xi }_2{\Xi }_4{\Xi }_9{e}^{-\lambda \tau }\\ +{\Xi }_2{\Xi }_4\lambda e^{-\lambda \tau }+{\Xi }_2{\Xi }_7{\Xi }_6{e}^{-\lambda \tau }+{\Xi }_4{\Xi }_8{\Xi }_3{e}^{-\lambda \tau }-{\Xi }_5{\Xi }_7{\Xi }_3{e}^{-\lambda \tau }+{\Xi }_7{\Xi }_3\lambda e^{-\lambda \tau }\end{array}$

$\begin{array}{l}\det (A-\lambda I)=[{\Xi }_1{\Xi }_5{\Xi }_9-{\Xi }_2{\Xi }_4{\Xi }_9-{\Xi }_1{\Xi }_8{\Xi }_6+{\Xi }_2{\Xi }_7{\Xi }_6+{\Xi }_4{\Xi }_8{\Xi }_3-{\Xi }_5{\Xi }_7{\Xi }_3\\ +({\Xi }_8{\Xi }_6+{\Xi }_2{\Xi }_4+{\Xi }_7{\Xi }_3-{\Xi }_5{\Xi }_9-{\Xi }_1({\Xi }_5+{\Xi }_9))\lambda +({\Xi }_5+{\Xi }_9+{\Xi }_1){\lambda }^2-{\lambda }^3]e^{-\lambda \tau }\text{ (114)}\end{array}$

We define for simplicity global parameters: ${\Omega }_0={\Xi }_1{\Xi }_5{\Xi }_9-{\Xi }_2{\Xi }_4{\Xi }_9-{\Xi }_1{\Xi }_8{\Xi }_6+{\Xi }_2{\Xi }_7{\Xi }_6+{\Xi }_4{\Xi }_8{\Xi }_3-{\Xi }_5{\Xi }_7{\Xi }_3$ ,

${\Omega }_1={\Xi }_8{\Xi }_6+{\Xi }_2{\Xi }_4+{\Xi }_7{\Xi }_3-{\Xi }_5{\Xi }_9-{\Xi }_1({\Xi }_5+{\Xi }_9)$ ; ${\Omega }_2={\Xi }_5+{\Xi }_9+{\Xi }_1$ ; ${\Omega }_3=-1$ .

Then we get $\det (A-\lambda I)=[{\displaystyle \sum _{k=0}^3{\Omega }_k}{\lambda }^k]e^{-\lambda \tau }$ ; $\det (A-\lambda I)=0$ , the general characteristics equation is $D(\lambda ,\tau )=[{\displaystyle \sum _{k=0}^{m=3}{\Omega }_k{\lambda }^k}]e^{-\lambda \tau }$ ; $Q_m(\lambda ,\tau )={\displaystyle \sum _{k=0}^{m=3}{\Omega }_k{\lambda }^k}$ . We study the occurrence of any possible stability switching resulting from the increase of value of the time delay $\tau$ for the general characteristics equation $D(\lambda ,\tau )$ . The characteristics equation is $D(\lambda ,\tau )=P_n(\lambda ,\tau )+Q_m(\lambda ,\tau )e^{-\lambda \tau }$ ; the expression for $P_n(\lambda ,\tau )$ is equal to zero: $P_n(\lambda ,\tau )=0$ . The expression for $Q_m(\lambda ,\tau )$ : $Q_m(\lambda ,\tau )={\displaystyle \sum }_{k=0}^m{\lambda }^k{q}_k(\tau )$ . The characteristics equation is $D(\lambda ,\tau )=P_n(\lambda ,\tau )+Q_m(\lambda ,\tau )e^{-\lambda \tau }$ where $P_n(\lambda ,\tau )=0$ , then ${\Omega }_k(\lambda ,\tau )=q_k(\tau )\forall k=0,1,2,3$ . We get characteristics equation is $D(\lambda ,\tau )=Q_m(\lambda ,\tau )e^{-\lambda \tau }$ . The $P_n(\lambda ,\tau )=0$ then the characteristic equation is reduced to $Q_m(\lambda ,\tau )e^{-\lambda \tau }=0$ . Since $e^{-\lambda \tau }\ne 0$ for any finite eigenvalue $\lambda$ therefore the characteristic equation can only be satisfied if $Q_m(\lambda ,\tau )=0$ . Thus, when $P_n(\lambda ,\tau )=0$ the delay term ( $\tau$ ) not determine the characteristic roots. The roots are simply the zeros of $Q_m(\lambda ,\tau )$ [24][25]. For any finite complex value of $\lambda$ , the exponential term $e^{-\lambda \tau }$ can never equal to zero ( $e^{-\lambda \tau }\ne 0$ ). Therefore, the roots of the system are entirely determined by the algebraic equation $Q_m(\lambda )$ . We can analyse stability using standard polynomial methods (like the Routh-Hurwitz criterion) on $Q_m(\lambda )=0$ . To analyse the stability of a third-degree characteristic polynomial $Q_m(\lambda ,\tau )={\displaystyle \sum }_{k=0}^{m=3}{\lambda }^k{q}_k(\tau )={\displaystyle \sum _{k=0}^{m=3}{\Omega }_k}{\lambda }^k$ , we use the Routh-Hurwitz criterion. A system is stable only if all roots

( $\lambda$ ) have negative real parts. The third-degree characteristic polynomial is $Q(\lambda )=\sum _{k=0}^{m=3}{\Omega }_k{\lambda }^k$ , $Q(\lambda )={\displaystyle \sum }_{k=0}^{m=3}{\Omega }_k{\lambda }^k={\Omega }_0+{\Omega }_1\lambda +{\Omega }_2{\lambda }^2+{\Omega }_3{\lambda }^3$ , we assume the leading coefficient ${\Omega }_3>0$ (since ${\Omega }_3=-1$ , we multiply the characteristic polynomial by (-1) and get the new global parameters: ${\Omega }_0=-{\Xi }_1{\Xi }_5{\Xi }_9+{\Xi }_2{\Xi }_4{\Xi }_9+{\Xi }_1{\Xi }_8{\Xi }_6-{\Xi }_2{\Xi }_7{\Xi }_6-{\Xi }_4{\Xi }_8{\Xi }_3+{\Xi }_5{\Xi }_7{\Xi }_3$ ,

${\Omega }_1=-{\Xi }_8{\Xi }_6-{\Xi }_2{\Xi }_4-{\Xi }_7{\Xi }_3+{\Xi }_5{\Xi }_9+{\Xi }_1({\Xi }_5+{\Xi }_9)$ ; ${\Omega }_2=-({\Xi }_5+{\Xi }_9+{\Xi }_1)$ ; ${\Omega }_3=1$ . Now leading coefficient ( ${\Omega }_3>0$ ), all other coefficients must be strictly positive: ${\Omega }_2>0$ , ${\Omega }_1>0$ , and ${\Omega }_0>0$ (necessary conditions). The sufficient condition (stable system) is ${\Omega }_2{\Omega }_1>{\Omega }_3{\Omega }_0$ [26][27][28]. If all conditions are fully met the system is stable. If any condition fails or a coefficient is negative ( ${\Omega }_k<0\forall k\in \left\{0,1,2,3\right\}$ ) then the system is unstable. We get marginally stable system if ${\Omega }_2{\Omega }_1={\Omega }_3{\Omega }_0$ while all coefficients are positive. A marginally stable system is defined as a dynamic system that neither grows infinitely nor settles down to zero when given a finite push or input. Instead, the system produces continuous, un-damped oscillations or a permanent, bounded offset that lasts forever. If we want to inspect the system stability from eigenvalues ( $\lambda \text{'}s$ ) then an unstable eigenvalue will have a positive or zero real part and a stable eigenvalue will have a negative real part. We inspect if any of the coefficients are negative or zero. If so, there is at least one value with a positive or zero real part which refers to an unstable node [29][30][31].