text:
\section{Introduction}

Kinetic Ising models have long been helpful in understanding general questions
of nonequilibrium statistical physics, as well as for elucidating particular issues
of experimental relevance \cite{Bray}. There is by now a vast body of research
studying these models under Glauber and Kawasaki dynamics \cite{Glauber,
Kawasaki}, the so called  models $A$ and $B$ respectively in the terminology of
Hohenberg and Halperin \cite{Hohenberg}. The first type of dynamic considers
single spin flip processes for describing relaxation towards equilibrium in a variety
of magnetic materials,  while the second case consists of spin exchanges specially
aimed to study stochastic processes under conservation of total magnetization.
This constraint makes relaxation slower,  and is instrumental in studying phase
separation, domain growth, and freezing observed after rapid cooling in several
systems, such as binary fluids and alloys \cite{Bray}. All these phenomena are
characterized by a long relaxation time $\tau$ which near critical points diverges as
$\tau \propto \xi ^z$, where $\xi$ denotes the spatial correlation length (proportional
to the typical system size), and $z$ is the dynamic exponent of the universality class
to which the dynamic belongs.

When it comes to one dimension ($1d$), these dynamics are still largely amenable
to experimental probe \cite{Privman} apart from being interesting in their own
right \cite{Godreche}, and as is well known, in the Glauber case allow for an exact
treatment \cite{Glauber, Privman}. In contrast to that latter aspect, the equations
of motion for the Kawasaki dynamics conform an infinite hierarchy which to our
knowledge do not admit analytic solution, so the problem can only be approached
from approximations \cite{Privman}. Even though there is no finite critical
temperature through which to quench, using random walk arguments to derive an
approximate kinetic equation for the kink density under ferromagnetic (F)
interactions \cite {Robin}, interestingly,  the dynamic exponent has been put
forward to be the same as observed in higher dimensions, namely  $z = 3$, after a
deep quench to a low temperature (LT). Although for $d \ge 2$ such exponent was
also accounted for by  renormalization group arguments \cite{Bray, Bray2} as well
as by surface dynamical considerations \cite{Huse}, other studies in $1d$ related
systems \cite{Luscombe} suggest instead a value of $z=5$ \cite{CST}, in
agreement with linear response schemes \cite{Zwerger}.

In this work we re examine this exponent by constructing and diagonalizing
numerically the kink evolution operator associated to the master equation of the
Kawasaki dynamics in finite chains. Following the physical picture given in
\cite{Robin}, at late evolution stages  and LT regimes such operator would
essentially describe a highly diluted system as the kink or domain wall density is
basically a measure of the average inverse domain size. Consequently, in evaluating
dynamic exponents from finite samples it might be expected that size effects will
pose no severe limitations for the kink representation.  In fact, as we shall see, the
kink density will get fairly small even for the low lying levels of the evolution
operator, that is what ultimately matters at large times. On the other hand there is
also the question about estimating dynamic exponents at finite temperatures for
which the  Kawasaki kinetics still remains critical,  as opposed to the exponential
decay of the Glauber dynamics in $1d$. Although in this situation dynamic
exponents are no longer related to domain growth - these are cut off by a finite
correlation length - they still provide the fundamental relation between the typical
system size and its relaxation time. We shall also address this issue by diagonalizing
the evolution operator of the original spins. Insofar temperature is not too low, these
latter become weakly correlated and finite-size effects will not be paramount.
A posteriori, our results will lend further support to this view.

Therefore these dual representations - kinks and spins - provide a means for probing
the robustness of dynamic exponents throughout different regimes. As a result, it
will turn out that for F couplings $z$ varies continuously from $\simeq 3.11$ to
$\simeq 2$ at low and high temperature (HT) scales respectively. Such
nonuniversality should come as no surprise, since in the HT limit the Kawasaki
dynamics simply reduces to a diffusive disordering kinetics. However under
antiferromagnetic (AF) exchanges these exponents are rather robust and, in line
with previous findings \cite{Luscombe}, their values will remain diffusive even at
LT regimes, just as the $z$ of the Glauber universality class. We shall return to this
question later on within a context of antikink operators.

The layout of this work is organized as follows. In Sec. II we recast the
$1d$-\,Kawasaki dynamics in terms of a quantum spin analogy that readily lends
itself to evaluate spectrum gaps of evolution operators, i.e. dynamic exponents,
both in the spin and domain wall representations. Due to detailed balance, either
of these descriptions can be brought to a symmetric representation by means of
simple non-unitary spin rotations. This simplifies considerably the subsequent
numerical analysis of Sec. III in which spectrum gaps are obtained via standard
recursive techniques \cite{Lanczos} in various situations. For now let us remark
that already modest chain lengths are able to yield clear finite-size trends both at
LT and HT scales. Finally,  Sec. IV contains a  summarizing discussion along with
some remarks on extensions of this work.

\section{Spin and kink representations}

As is well known, the Ising model has no intrinsic dynamics because all spin
operators involved in its Hamiltonian commute with one another. Therefore, an
{\it ad hoc} dynamics must be prescribed by coupling the system to a heath bath
at temperature $T$ so as to induce energy changes in the model. This is described
by a master equation \cite{Kampen} for which we introduce briefly here some
preliminary considerations.

Basically, the dynamics is associated to a gain-loss equation, constructed
generically as
\begin{equation}
\label{master}
\partial_t \,P(s,t) = \sum_{s'} \left[\: W (s'\to s)\, P(s',t)\,
- \, W (s \to s') \,P(s,t)\:\right]\,,
\end{equation}
which governs the evolution of the probability $P (s,t)$ that the system will be at
state $\vert s \rangle$ at time $t$. The elementary change steps are embodied in
the transition probability rates $W ( s \to s')$ per unit time at which configuration
$\vert s \rangle\,$ evolves  to $\vert s' \rangle\,$. For our purposes, it is convenient
to think of  this equation as a Schr\"odinger evolution in an imaginary time, namely
$\partial_t \vert P (t)\, \rangle = - H \vert P (t) \,\rangle$, under a pseudo Hamiltonian
or evolution operator $H = H_d + H_{nd}$ whose diagonal and non-diagonal matrix
elements are given by
\begin{eqnarray}
\label{diag}
\langle\,s\,\vert\, H_d\,\vert\,s\,\rangle &=& \sum_{s'\ne s}\, W(s \to s')\,,\\
\label{non-diag}
\langle\,s'\,\vert\,H_{nd}\,\vert\,s\,\rangle &=& -\,W(s \to s')\,.
\end{eqnarray}
Formally, this enables one to derive all subsequent probability distributions
$\vert P(t) \,\rangle \equiv \sum_s P (s,t)\, \vert s \rangle\,$ from the action of $H$
on a given initial distribution, i.e. $\vert P(t) \,\rangle = e^{- H\,t} \vert P(0)\,
\rangle$ \cite{Kampen}. In particular, the relaxation time of any observable is
singled out by the eigenvalue $\lambda_1$ corresponding to the first excitation
mode of $H$, i.e. $\tau^{-1} = {\rm Re}\,\lambda_1 > 0$, whereas the steady
state merely corresponds to an eigenvalue $\lambda_0 = 0$\, \cite{Kampen}.
If the steady configuration should actually coincide with the Boltzmann equilibrium
distribution, the above rates must be constrained by detailed balance, that is
$W (s \to s') \, e^{- \beta E (s)} =  W (s' \to s) \,e^{- \beta E (s')}$,  $\forall \,s,\,s',$
where $E$ stands for the respective energy configurations of the system in question,
and $\beta \equiv 1/k_B T$.

Detailed balance itself can not determine entirely the form of such rates, thus for
the specific case of the Kawasaki dynamics hereafter considered we take up the
common choice \cite{Bray,Kawasaki}
\begin{equation}
\label{rates0}
W (s \to s') = \frac{1}{2}\Big[\,1- \tanh \Big(\,\frac{\beta }{2} \,
\Delta E_{s, s'} \Big)\,\Big]\,,
\end{equation}
where $\vert s \rangle ,\, \vert s'\rangle $ are states differing at most in a pair of
nearest neighbor (NN) spins exchanged at some location $i = 1,2,\,...\,L$, and
$\Delta E _{s , s'} \equiv E (s') -  E (s)$ is the change of the Ising energy $E = - J
\sum_i s_i \,s_{i+1},\, (s_i = \pm 1)$ either with F $(J > 0)$, or AF $(J < 0)$
interactions. Therefore depending on the spin states at locations $i-1$ and $i+2$,
the rate at which spins $(s_i,\, s_{i+1}) = (s,\,-s)$ exchange their states results in
\begin{equation}
\label{rates}
R_{i-1,i+2} (\pm K) = \frac{1}{2} \pm \frac{\alpha_{_K}}{4} (\,s_{i+2} - s_{i-1})
\,,
\end{equation}
where $K = \beta J$, $\alpha_{_K} = \tanh 2 K$, and the signs $\pm$ denote
forward and backward hoppings as depicted schematically in Fig.\,1.
\begin{figure}[htbp]
\vskip -8.5cm
\centering
\includegraphics[width=0.8\textwidth]{F1.ps}
\vskip -9cm
\caption{Transition probability rates for neighboring spin
exchanges, as defined in Eq.\,(\ref{rates}).}
\end{figure}
If we think of the spin configurations $\vert s \rangle$ as being already diagonal
in the $z$-direction, say, then by promoting $s_i$ variables to Pauli matrices
$\sigma_i^z$, the operational analog of Eq.\,(\ref{non-diag}) will read
\begin{equation}
\label{nondiag-spins}
H_{nd} = - \sum_i  \left[\, R^z_{i-1,i+2} (K) \,\sigma^+_{i+1}\sigma^-_i + \,
R^z_{i-1,i+2} (-K)\, \sigma^+_i \sigma^-_{i+1}\right] \,,
\end{equation}
where $\sigma^+, \sigma^-$ are the usual spin-$\frac{1}{2}$ raising and
lowering operators, while $R^z_{i-1,i+2} (\pm K) = \frac{1}{2} \pm
\frac{\alpha_{_K}}{4} (\,\sigma^z_{i+2} - \sigma^z_{i-1}) \,,$ is simply the
operational counterpart of Eq.\,(\ref{rates}).  As for the diagonal elements of
(\ref{diag}) needed for conservation of  probability, notice that these basically count
the number of hoppings in which a given configuration can evolve to different ones
by exchanging NN spins at a time.  This can be properly tracked down in terms of
number operators $\hat n = \sigma^+ \sigma^- $ probing and weighting all NN
spins, namely
\begin{equation}
\label{diag-spins}
H_d = \sum_i \left[\,  R^z_{i-1,i+2} (K) \,\hat n_i\, (1- \hat n_{i+1})
+ R^z_{i-1,i+2} (-K)\,\hat n_{i+1}\, (1- \hat n_i)\,\right]\,.
\end{equation}

Although the correlated exchange terms of Eq.\,(\ref{nondiag-spins}) leave us with
a non-symmetric evolution operator, in preparation for the numerical analysis of
Sec. III  we can make some progress by exploiting detailed balance. This latter
warrants the existence of representations in which $H$ is symmetric and thereby
diagonalizable \cite{Kampen}. For our purposes, it suffices to consider the diagonal
non-unitary similarity transformation
\begin{equation}
\label{trans1}
S = e^{-\frac{K}{2}\,\sum_j \sigma^z_j \, \sigma^z_{j+1} }\,,
\end{equation}
under which the hopping terms of $H_{nd}$ transform as
\begin{equation}
\label{transform}
\sigma^{\pm}_{i+1} \sigma^{\mp}_i \to \exp \left[\,\mp K \left( \sigma^z_{i+2}
- \sigma^z_{i-1 }\right)\, \right]\, \sigma^{\pm}_{i+1} \sigma^{\mp}_i\,,
\end{equation}
while leaving invariant all the above mentioned diagonal operators. Hence, after
straightforward manipulations it can be readily verified that ${\rm H}_{spin} =
S  H S^{-1}$ actually produces a self adjoint spin representation which for periodic
boundary conditions (PBC) is found to be
\begin{eqnarray}
\nonumber
{\rm H}_{spin} &=& - \frac{1}{8} (\,1 + \gamma_{_K} ) \, \sum_i \, \left(\, 1 +
\tanh^2 \! K \,\sigma^z_{i-1}\sigma^z_{i+2}\right)
 \left( \, \sigma^x_i \sigma^x_{i+1} + \sigma^y_i \sigma^y_{i+1}\,  \right)\\
\label{hermitian-spin}
&+& \frac{1}{4} \sum_i \left[\,1 \,+\,\alpha_{ _K} \,\sigma^z_i \sigma^z_{i+2}
 \,-\, (1+ \alpha_{ _K}) \,\sigma^z_i \sigma^z_{i+1} \, \right]\,,
\end{eqnarray}
where $\gamma_{_K} = {\rm sech}\, 2 K$.
In passing, it is worth pointing out that in the HT limit this operator coalesces into a
Heisenberg ferromagnet, i.e. ${\rm H}_{spin} \to -\frac{1}{4} \sum_i  \left(\,\vec
\sigma_i \cdot \vec\sigma_{i+1} - 1  \right)$, with a spectrum gap behaving like
$\propto 1/L^2$, as it should for a disordering diffusive kinetics. On the other hand,
in such regime the equilibrium correlation length becomes of the order of the lattice
spacing, thus the numerical diagonalization of ${\rm H}_{spin}$ in finite chains
can be a reasonable starting point to evaluate dynamic exponents as temperature
is quenched slightly (though {\it instantly} to a constant value, because transition
rates have been taken time independent  throughout).

When it comes to LT regimes however, Eq.\,(\ref{hermitian-spin}) is not of much
practical use, since the  system approaches a state consisting of long magnetic
domains whose typical sizes are of the order of the correlation length $\xi \sim
e^{2 \vert K \vert}$ (nonetheless, for $J > 0$ see the upper bounds to $z$ and their
extrapolations provided in Sec. III A). In this more interesting scenario the dynamics
is strongly affected by whether the coupling exchanges are F or AF, thus it is
convenient at this point to separate the discussion accordingly.

\subsection{Kink representation: $J > 0$}

As was mentioned in Ref.\,\cite{Robin}, the dynamics of this situation is basically
mediated by single spins detaching from a domain wall (a kink),  and then
performing a random walk. Some of  these spins will return to their original
domains, while others will reach the next ones. On average, as a result of many
of these processes the domains themselves will perform a random walk until they
meet and merge into larger domains. Besides, single meandering spins may
eventually nucleate and trigger the growth of a new domain. The key issue in this
description is that for large times very few kinks will survive at LT scales.
In particular, in equilibrium these are totally uncorrelated and with an average
density $\sim 1/\xi$, whereas that of the spins mediating the whole process is
of order $1/\xi^2$ as they are made up of two consecutive kinks (see
Fig.\,\ref{Fkinks}a below).  Thus, for large times it makes sense to attempt a
nonequilibrium description in terms of kinks rather than spins. Presumably
(as confirmed later on in Sec. III), the first excitation mode of the kink evolution
operator will also correspond to a diluted state, so size effects in a numerical
estimation of $z$ will not be as severe as in the spin description.

\begin{figure}[htbp]
\vskip -6.75cm
\centering
 \includegraphics[width=0.65\textwidth]{F2.ps}
\vskip -5.25cm
\caption{Transition rates for (a) diffusion of double kinks, and (b) assisted
deposition-evaporation of kinks. These latter are denoted by vertical lines
separating domains of opposite spin orientations.}
\label{Fkinks}
\end{figure}

To recover the equilibrium behavior of the Ising model from now on we consider
the case of zero magnetization, so this dual representation corresponds to
the two-to-one mapping outlined schematically in Fig.\,\ref{Fkinks}. If we imagine
these kinks as hard core particles $A$, then the Kawasaki dynamics involve two
basic processes (i) dimer diffusion $A+ A + \emptyset \rightleftharpoons
\emptyset + A + A$, and (ii) `assisted' deposition-evaporation $\emptyset + A +
\emptyset \rightleftharpoons A + A + A$,  as schematized respectively in
Figs.\,\ref{Fkinks}a and \ref{Fkinks}b. The first situation represents the
meandering spins referred to above and involves no energy changes (rates $1/2$).
Notice that these dimers (double kinks) do not preserve their identity as they may
eventually dissociate when contacting a domain wall (single kink). The second
process corresponds to a spin detachment (attachment) from a domain wall, and
its rates $\epsilon_{_ K}$  ($\epsilon_{_ {\! - K}}$) parallel those of
Eq.\,(\ref{rates}) when $s_{i-1} = - s_{i+2}$, thus
\begin{equation}
\epsilon_{_ {\!\pm K}} = \frac{1}{2}\,\left(\,1 \mp \tanh 2K\, \right)\,.
\end{equation}

Therefore, following the reasoning steps discussed for the spin representation
and reinterpreting kinks as $\frac{1}{2}$-spinors in a dual chain, clearly the
non-diagonal terms  $H^a_{nd} + H^b_{nd}$ of the evolution operator associated
to the (a) and (b) processes of Fig.\,\ref{Fkinks} can be constructed as
\begin{eqnarray}
\label{double-kink}
H^a_{nd} &=& -\frac{1}{2} \sum_i  \hat n_i \left(\, \sigma^+_{i+1}
 \sigma^-_{i-1} + {\rm h.c.} \right)\,,\\
\label{dep-eva}
H^b_{nd} &=& -\sum_i \hat n_i \left(\,\epsilon_{_K}\, \sigma^+_{i-1}
 \sigma^+_{i+1} + \, \epsilon_{_{\!-K}}\, \sigma^-_{i-1} \sigma^-_{i+1} \right)\,,
\end{eqnarray}
with number operators acting here as projectors that rule out vacancy mediated
processes. To recast the non-diagonal operator into a symmetric representation we
recur once more to detailed balance and rotate all $\sigma^{\pm}$\,'s around the
$z$  direction using a  common pure imaginary angle $\varphi = i K$.  This is
carried out by the non-unitary similarity transformation
\begin{equation}
\label{i-rotation}
{\rm S} = e^{ \frac{K}{2}\,\sum_j \sigma^z_j}\,.
\end{equation}
for which it is straightforward to show that
\begin{equation}
\sigma^{\pm}_{i-1} \sigma^{\pm}_{i+1}  \to
e^{\pm 2 K} \sigma^{\pm}_{i-1} \sigma^{\pm}_{i+1}\,,
\end{equation}
thus producing the symmetrization of  $H^b_{nd}$ while keeping  all terms of
$H^a_{nd}$ (already hermitian) unaltered.

As before, preservation of probability is taken into account by a diagonal operator
$H^a_d + H^b_d$ balancing both of the above (mutually exclusive) events, in turn
probed respectively by
\begin{eqnarray}
H^a_d &=& \frac{1}{2} \sum_i \left[\,\hat n_{i-1}\, \hat n_i  \left( 1 - \hat n_{i+1}
\right) \, +\,  \left( 1 - \hat n_{i-1}\right) \hat n_i \, \hat n_{i+1}\, \right]\,,\\
H^b_d &=& \sum_i \left[\,\epsilon_{_K} \left(1-\hat n_{i-1}\right)\, \hat n_i \,
\left(1-\hat n_{i+1}\right)\,  + \, \epsilon_{_{\!-K}} \, \hat n_{i-1}\, \hat n_i \,
\hat n_{i+1} \, \right]\,.
\end{eqnarray}
Since the imaginary rotation (\ref{i-rotation}) is also diagonal, it has no effect
on these latter operators. Thus, collecting all terms and using PBC throughout,
after some algebraic steps we finally obtain a self adjoint representation
for the kink evolution operator, namely
\begin{eqnarray}
\nonumber
{\rm H}_{kink} &=&  - \, \frac{1}{8} \sum_i \, \big(\,1+\sigma^z_i \,\big)\,
\left[\, \left (\,1+ \gamma_{_K}\right) \,
\sigma^x_{i-1}\, \sigma^x_{i+1}  + \left (\,1- \gamma_{_K}\right) \,
\sigma^y_{i-1}\, \sigma^y_{i+1}\, \right] \\
\label{hermitian-kink}
&+& \frac{1}{4} \sum_i  \left(\,1 + \alpha_{_K} \, \sigma^z_i\,  \sigma^z_{i+1}
 \right) + \frac{1}{2} \,\epsilon_{_{-K}} \sum_i \sigma^z_i \;\; , \;\; K > 0\,.
\end{eqnarray}
As can be readily verified from this equation, in the limit $T \to 0^+$ the action
of ${\rm H}_{kink}$ on any configuration having non NN kinks vanishes like
$O(e^{-4 K})$, thus yielding a metastable state. In the spin representation this
corresponds to configurations with domains lengths larger than the lattice spacing,
which amounts to hindering spins to detach and so diffuse through. In turn, from
Eq.\,(\ref{hermitian-spin}) it can be checked that this situation also yields
metastable states of ${\rm H}_{spin}$ so long as $J > 0$. Clearly, the number of
these configurations grows exponentially with the  system size, and most of them
take over the asymptotic dynamics with long lifetimes $\propto \xi^2$. This marks
an important difference with respect to the AF dynamics to be introduced briefly
in Sec. II B.

As an aside, we finally mention that by construction ${\rm H}_{kink}$ not only
preserves the parity of kinks $e^{i \pi \sum_j \hat n_j}$ (being even for PBC),
but in turn satisfies
\begin{equation}
\Big[\, {\rm H}_{kink}\,,\, \sum_j \sigma^z_j \,e^{i \pi \sum_{k \le j}\hat n_k }\,
\Big] = 0 \,,
\end{equation}
which simply expresses the conservation of the total spin magnetization in the
original system. In practice, for the numerical evaluation of gaps we will just  build
up the adequate kink states from the corresponding spin ones.

\subsection{Antikinks: $J < 0$}
While most of the above ideas and procedures applies to the AF dynamics as well,
ultimately the relaxation to equilibrium becomes faster than in the F case.  By
changing from a description based on kinks to AF domain walls, i.e. to antikinks, we
can easily construct an evolution operator that for large times and LT
regimes essentially describes a diluted antikink system.  As before, if we think of
antikinks as hard core particles $A$ (former vacancies $\emptyset$ under F
couplings), the basic processes now involve (i) next NN hopping, i.e. $A + \emptyset
+ \emptyset \rightleftharpoons  \emptyset + \emptyset + A$, and (ii) deposition -
evaporation `assisted'  by vacancies, that is $\emptyset + \emptyset + \emptyset
\rightleftharpoons A + \emptyset + A$.

\begin{figure}[htbp]
\vskip -8cm
\centering
 \includegraphics[width=0.65\textwidth]{F3.ps}
\vskip -4cm
\caption{Transition rates for (a) next nearest neighbor antikink hoppings mediated
by vacancies, and (b) assisted  deposition-evaporation of antikinks. As in Fig. 2,
these latter are represented by vertical lines but separating antiferro domains.}
\label{AFkinks}
\end{figure}

With the aid of these schematic events, illustrated respectively in
Figs.\,\ref{AFkinks}a and \ref{AFkinks}b, and after carrying out the rotation
referred to in Eq.\,(\ref{i-rotation}) but with an argument $\bar\varphi = -i K$,
the symmetric representation of the antikink operator can be finally cast as
\begin{eqnarray}
\nonumber
{\rm H}_{anti} &=&   \frac{1}{8} \sum_i \,  \big(\,1-\sigma^z_i \,\big)\,
\left[\, \left (\,1+ \gamma_{_K}\right) \,
\sigma^x_{i-1}\, \sigma^x_{i+1}  + \left (\,1- \gamma_{_K}\right) \,
\sigma^y_{i-1}\, \sigma^y_{i+1}\, \right]  \\
\label{hermitian-antikink}
&+& \frac{1}{4} \sum_i  \left(\,1 + \alpha_{_K} \, \sigma^z_i \,
\sigma^z_{i+1} \right) - \frac{1}{2} \,\epsilon_{_{-K}} \sum_i \sigma^z_i
\;\; , \;\; K < 0\,.
\end{eqnarray}
Of course, by exchanging the role of particles and vacancies, i.e. using the canonical
transformation $\sigma^{\pm} \to \sigma^{\mp}$, one recovers the form of
Eq.\,(\ref{hermitian-kink}), but this likelihood is only superficial. Note that for
$K \to -\, \infty$ the effect of the uniform field $\epsilon_{_{-K}}$ vanishes like
$O (e^{4K})$, so the {\it only} jammed state on which the action of the diagonal
terms cancels out at this order is the vacuum of antikinks. Thus, in the limit
$T \to 0^+$ the dynamics is no longer dominated by metastable states (instead
proliferating exponentially for $J > 0$). In this situation AF domains can grow
unhindered, as single antikinks now detach and diffuse freely (Fig.\,\ref{AFkinks}),
thus allowing the dynamics to run smoothly towards equilibrium. As a result the
dynamic exponent, whose numerical analysis we next turn to consider, becomes
smaller (actually diffusive) than for F exchanges.

\section{Numerical Results}

We now investigate numerically the spectral gaps of the Kawasaki operators
given in Eqs.\,(\ref{hermitian-spin}), (\ref{hermitian-kink}), and
(\ref{hermitian-antikink}) along with their dynamic exponents in different
situations. As a consistency check  first we verified that the (rotated) Boltzmann
distribution corresponds in fact to the steady state of those operators. This also
served to start up the Lanczos recursion with a random state but chosen orthogonal
to that equilibrium configuration.  Thereafter we obtained the first excited
eigenmodes of our symmetric ${\rm H}'s$ using periodic chains of up to 24 sites,
the main limitation for this being the exponential growth of the space
dimensionalities.  Another restrictive issue is that below  $k_B T / \vert J \vert \sim$
0.1 - 0.2\, the Lanczos convergence slows down progressively because in most
situations the spacing of low lying levels turns out to decay as $e^{-4 \vert K
\vert }$. Thus, in what follows we content ourselves with giving results above that
region where nonetheless clear saturated tendencies can be already obtained.

\subsection{$J > 0$}

Turning to the F dynamics, in the insets of Figs.\,\ref{F-exponents}a and
\ref{F-exponents}b we exhibit the finite-size behavior of spectral gaps for the case
of spins and kinks respectively. They are all consistent with a gap vanishing  like
$1/L^z$, although with a nonuniversal temperature dependent dynamic exponent
$z\,(T)$.  For the purpose of observing in more detail the trend of size effects on
these exponents under different temperature regimes, we considered a sequence
of effective approximants defined as
\begin{equation}
\label{approximants}
z_L = \frac{\ln\left[ \,\lambda_1 (L-2) / \lambda_1 (L) \,\right]}
{\ln \left[ \,L / (L-2) \,\right] }\,,
\end{equation}
($L$ even), which simply provide successive measures of the gap closing in either
of the above representations. These are shown in the main panels of
Figs.\,\ref{F-exponents}a and \ref{F-exponents}b. As expected, in the HT
region ${\rm H}_{spin}$ is able to yield fairly convergent estimations
of $z\, (\,\sim 1.99)$, so recapturing the plain diffusive limit above  $T / J \simeq 3$
(henceforth the Boltzmann constant $k_B$ is set equal to 1). Besides, in this
representation the above approximants turn out to constitute upper bounds
of this exponent for almost all temperatures studied, though exhibiting
different spreadings. We attempted to extrapolate these bounds using several
procedures, but owing to the variable spreading this resulted in a rather noisy
limiting curve. Nevertheless, in the region that most interests us the spreading
becomes more stable, thus upon using a van den Broeck - Schwartz extrapolation
scheme \cite{VBS} we found that  $z \sim \, 3.1(0)$, in reasonable agreement with
analysis and Monte Carlo simulations under instantaneous quenches \cite{Robin}.
Moreover, this result is also consistent with estimations of $z$ arising from
${\rm H}_{kink}$ approximants around the same region ($T / J \alt 0.4$) where,
as discussed above, the kink description in finite sizes is most reliable. In fact, no
extrapolations are needed here since already the values of $z_{20} \alt z_{22}
\alt z_{24}$ follow very closely one another, collapsing near a value of $z \simeq
3.11$. This is the main result of this section.

As an added bonus, note that in the LT regime the kink approximants
come out conforming a sequence of lower bounds which conveniently complement
the upper ones provided by the spin representation. This is illustrated in
Fig.\,\ref{F-exponents}c where for $ T / J \alt 0.4$ both kink and spin approximants
should actually enclose the value of $z$ in the limit $L \to \infty$ (presumably much
closer to the kink border). Using the results of our maximum reachable sizes, this
suggests that in the ordering regime the dynamic exponent is bounded as $3.11
\alt  z < 3.18$. For $0.5 \alt T / J \alt 2$ the trend of finite sizes of the kink
approximants now reverses and, alike the spin  $z_L$'s, turn out to yield upper
bounds of $z$. The kink approximants converge slightly faster in this region (see
inset of Fig.\,\ref{F-exponents}c), until for $T/J \agt 2$ the trend of bounds reverses
once more. Hence, at HT scales the actual value of $z$ is enclosed again by our kink
and spin approximants though this time these latter take the lead, as it should.

\begin{figure}[htbp]
\vskip -3cm
\hskip -2.5cm
\includegraphics[width=0.65\textwidth]{F4a.ps}
\hskip -3cm
\includegraphics[width=0.65\textwidth]{F4b.ps}
\vskip -6cm
\hskip -2cm
\includegraphics[width=0.65\textwidth]{F4c.ps}
\vskip -3.5cm
\caption{(Color online) Effective dynamic exponents of Eq.\,(\ref{approximants}) for
$J > 0$ using $12 \le L \le 24$ with $L$ even, computed for (a) spins, and (b) kinks.
Sizes increase  downwards almost  throughout (a), and upwards in (b) at both
leftmost and rightmost  temperature regions. The respective insets display the gap
behavior with the system size for $T/ J = 10, 1, 1/2$ and  $1/4$ (top to bottom).
To compare slopes, data have been shifted upwards with respect to $T/ J = 10$.
Straight lines are fitted with slopes  $-z_{24}(T)$ as calculated from the main
panels. In (c) we compare these exponents for kinks and spins (dashed and solid
lines). At low as well as at high temperatures $z_k$ and $z_s$ yield respectively
lower and upper bounds for the thermodynamic limit of $z$. The inset of (c)
provides a measure of convergence for these quantities by depicting $\Delta_L z
\equiv z_{L-2} - z_L $ for our maximum sizes. Horizontal short doted line in (a)
stands for the extrapolations referred to in the text.}
\label{F-exponents}
\end{figure}

As mentioned earlier, the LT convergence of kink approximants was somehow
foreseeable on the basis of equilibrium considerations but it is not yet clear whether
these apply to nonequilibrium as well. In Fig.\,\ref{kink-density} however we show
that this also the case, at least for the sizes at hand. There, we display the kink
density $\rho_1 = \frac{1}{L} \sum_i \langle \psi_1 \vert \hat n_i \vert \psi_1
\rangle$ for the first excited mode of ${\rm H}_{kink}$ which for $T/J \alt 1$ closely
follows the diluted kink picture already expected for equilibrium, i.e. $\rho_1 \to 2/L$.
Also, preliminary evaluations of kink-kink correlations $\langle \psi_1 \vert \hat n_i
\hat n_j \vert \psi_1 \rangle$ show that these become negligible below $T/J \sim 0.5$,
thus resembling the fully segregated equilibrium state. In passing it is worth
pointing out that while $\lambda_1$ is doubly degenerate, the correlations and
$\rho_1$ remain indistinguishable in both of these eigenmodes.
\begin{figure}[htbp]
\centering
\vskip -0.5cm
\includegraphics[width=0.472\textwidth]{F5.ps}
\caption{Density of kinks for  the first excitation (solid line), and ground state
(dashed line) of ${\rm H}_{kink}$ [\,Eq.\,(\ref{hermitian-kink})\,] with $L = 24$ and
$J > 0$. Lowermost dots indicate the equilibrium density $(2/L$) when $T \to 0$.}
\label{kink-density}
\end{figure}

\vspace{-0.5cm}

\subsection{$J < 0$}

In studying numerically the AF dynamics, the approximants (\ref{approximants})
pose severe difficulties at intermediate temperatures. Owing to low lying level
crossings, now the spectrum gap is bisected in temperature regions having opposite
monotonic behaviors. For instance, it is instructive to consider the case of $L=4$
where this feature already appears in the spectrum of ${\rm H}_{spin}$. By
diagonalizing its associated $6 \times 6$ stochastic matrix, we find that there
are two branches of eigenvalues crossing at $K^* = - \frac{1}{4} \ln 2$, thus
producing the non-analytic gap
\begin{equation}
\label{L=4}
\lambda_1 (4) = \cases{2\, \epsilon_K\;\;\,, \;{\rm threefold \;\;
                                   degenerate\;\; if}\;\; K \ge K^*\,,\cr
                                    4 \, \epsilon_{-K}\,, \; {\rm nondegenerate\;\;
                                    otherwise}\,,}
\end{equation}
as seen in the uppermost curve of Fig.\,\ref{lambdas}a. Alike this simple case, for
larger sizes it is found numerically that  $\lambda_1 (L)$ remains degenerate only
above certain temperatures below which level crossings occur, and $\lambda_1$
results in a nondegenerate value. This is signaled by the emergence of the cusps
observed in Fig.\,\ref{lambdas}a.  At LT regions the gap recovers the exponential
decay referred to above but the data collapse now precludes finite-size estimations
of effective exponents from Eq.\,(\ref{approximants}). Although level crossings
come out at successively smaller temperatures, these can not keep pace with the
increasing lattice sizes. To bypass these limitations we resorted to the antikink
operator especially constructed for this region [\,Eq.\,(\ref{hermitian-antikink})\,].
As can be seen in Fig.\,\ref{lambdas}b, now the LT spectrum no longer collapses so
the approximants (\ref{approximants}) can be employed once more. However, as
evidenced by the cusps of that figure, the problem of low lying level crossings yet
persists. In nearing the crossing temperatures this brings about pronounced size
effects in $z_L$, partly because of the mismatches occurring between eigenvalue
branches, alike those observed in Fig.\,\ref{lambdas}a.

\begin{figure}[htbp]
\hskip -0.5cm
\includegraphics[width=0.472\textwidth]{F6a.ps}
\hskip 0.7cm
\includegraphics[width=0.472\textwidth]{F6b.ps}
\vskip 0.5cm
\caption{Spectrum gaps of (a) spin, and (b) antikink operators
[\,Eqs.\,(\ref{hermitian-spin}) and (\ref{hermitian-antikink}) respectively\,], for
various inverse temperatures $K = J/T$, ($J < 0$).  Sizes increase from top to
bottom. In (a) the dashed line stands for the soluble case $L = 4$ already indicating
the appearance of a cusp due to level crossings [\,see Eq.\,(\ref{L=4})\,], a feature
holding also for $L = 6,8,...,24$. At low temperatures $\lambda_1 (L)$ decays as
$e^{- 4\vert K \vert}$ with an almost  size independent amplitude. In (b) this latter
problem is remedied although level crossings for $8,10, ... \,,$ and 24 sites now give
rise to non-analiticities near which size effects become severe
(cf. Fig.\,\ref{a-exponents}).}
\label{lambdas}
\end{figure}

All the above considerations result in fair convergent estimations of dynamic
exponents both for LT and HT regimes, but as displayed in
Fig.\,\ref{a-exponents}, the intermediate zone is out of reach. Despite this setback,
what matters is the usefulness of the antikink representation
(\ref{hermitian-antikink}) in the more interesting LT region which
otherwise would have remained inaccessible. As in the F case, no extrapolations to
the thermodynamic limit are needed here since already for $T/\vert J \vert \alt 0.3$
our higher antikink approximants lie within the interval $1.98 \simeq z_{20} \alt
z_{22} \alt z_{24} \simeq 1.99$, thus indicating a Glauber universality class for
this ordering regime, cf. \cite{Luscombe}.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.472\textwidth]{F7.ps}
\vskip 0.5cm
\caption{Effective dynamic exponents as defined in Eq.\,(\ref{approximants}),
for Hamiltonians (\ref{hermitian-spin}) and (\ref{hermitian-antikink})
with $J < 0$ and $L = 10, 12,\,...\,, 24$ increasing upwards. In nearing the
cusps of Fig.\,\ref{lambdas} size effects become progressively pronounced,
i.\,e. $\lambda_1 (L)$ and $\lambda_1 (L-2)$ belong to different branches.
High and low temperature regimes are both consistent with a typical
diffusive exponent $( z_{24} \simeq 1.99)$. }
\label{a-exponents}
\end{figure}
In the other extreme, when disorder prevails, already for $T/\vert J \vert \agt 3$
the spin approximants derived from the data of Fig.\,\ref{lambdas}a converge
within one percent to our LT exponents, thus suggesting for these
latter a rather robust diffusive behavior.

Finally, in Fig.\,\ref{antikink-density} we exhibit the antikink density $\rho_1$ for
the first excited state of ${\rm H}_{anti}$. As before, this confirms the diluted
antikink picture expected at LT scales, in turn accounting for the convergent
approximants obtained below $T/J  \alt 0.3$. However, this time $\rho_1$ does not
quite follow there the vanishing equilibrium density of the plain AF vacuum state
and adopts instead the limiting value referred to in Fig.\,\ref{kink-density},
indicating the presence of two antikink excitations.  At intermediate temperatures
the jump of $\rho_1$ just reflects the aforementioned level crossings.

\begin{figure}[htbp]
\centering
\vskip -0.1cm
\includegraphics[width=0.472\textwidth]{F8.ps}
\caption{Density of antikinks for the first level and ground state of  operator
(\ref{hermitian-antikink}) with $L = 24$ and $J < 0$ (solid and dashed lines
respectively). For comparison, horizontal dots indicate the equilibrium  density
of ferro kinks $(J > 0)$ at $T=0$, cf. Fig.\,\ref{kink-density}. The discontinuity
denoted by vertical dots stems from level crossings in the low lying spectrum of
${\rm H}_{anti}$ (also, see cusp in the lowermost curve of Fig.\,\ref{lambdas}b).}
\label{antikink-density}
\end{figure}

\section{Concluding discussion}

To summarize,  we have constructed symmetric representations of the Kawasaki
dynamics in the Ising chain using a quantum spin analogy with both direct an dual
processes. The kink description \cite{Robin}, either in its F or AF versions
[\,Eqs.\,(\ref{hermitian-kink}) and (\ref{hermitian-antikink})\,], is well suited to low
temperature regimes where it is able to provide a rather fast finite-size convergence
to dynamic exponents (Figs.\,\ref{F-exponents}b and \ref{a-exponents}). For F
couplings, the ordering kinetics arising from these latter ($z \simeq 3.11$) tuns out
to be slightly slower than the Lifshitz - Slyozov behavior [\,$\xi (t) \propto t^{1/3}$]
characteristic  of higher dimensions \cite{Huse}. By contrast, for AF interactions
the kinetics is no longer activated by metastable states, and dynamic exponents
converge rapidly ($z \simeq 1.99$) to the typical diffusive values of the Glauber
universality class [\,$\xi (t) \propto t^{1/2}$].

On the other hand, the spin representation [\,Eq.\,(\ref{hermitian-spin})\,] is more
appropriate for higher temperature scales where the kinetics remains slow and is
still characterized by relaxation times diverging as $L^z$ (inset of
Fig.\,\ref{F-exponents}a). For the F dynamics, the finite-size approximants of these
exponents turn out to form a sequence of upper bounds which also extends down to
low temperature regions (Fig.\,\ref{F-exponents}a), where it nicely complements the
series of lower bounds emerging from the kink approach (Fig.\,\ref{F-exponents}b).
On wider temperature scales, both representations are consistent with a
nonuniversal set of dynamic exponents interpolating continuously between a
subdiffusive ordering kinetics and the plain diffusive limit (Fig.\,\ref{F-exponents}c).
For $J < 0$, part of this intermediate region proved inaccessible due to size effects
caused by low lying level crossings in the spectrum of both spin and antikink
operators (cusps of Fig.\,\ref{lambdas}). However, as the dynamics of this case is
not activated, the diffusive exponents obtained by both representations within their
natural range of applicability (Fig.\,\ref{a-exponents}) suggest a rather universal
(Glauber) dynamics throughout.

The success of the kink description at low temperature regimes would probably
allow to extend the ideas of this work to other interesting $1d$-dynamics under
instantaneous quenches (i.e. constant transition rates). For instance, the
nonuniversal behavior alleged for the ordering kinetics of the alternating bond
Kawasaki chain \cite{Luscombe}, might well be further investigated with our
numerical approach.  Other lines of research that would also be worth pursuing are
certain dynamical processes involving composite objects (e.g. dimers) which exhibit
strongly broken ergodicity as a result of having an extensive number of conservation
laws \cite{Barma}. As the issue of universality classes in nonequilibrium statistical
systems is often linked to the existence of these latter \cite{Odor}, it would be
important to determine whether dynamic exponents actually depend on the
subspaces where the evolution takes place.  Dual representations constructed
appropriately for such processes could therefore provide reliable computations
of these exponents within ordering regimes. Further work in this direction is in
progress.

\section*{Acknowledgments}

The author wishes to thank M. Arlego, G. L. Rossini, and F. A. Schaposnik for
helpful comments and discussions. Support of CONICET, Argentina under Grants
No. PIP 1691 and No. PICT ANCYPT 20350, is acknowledged.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

The quest to find solutions to energy storage for the renewable industry and the rise of mobile robots, electric vehicles, and drones have led to further development of new generations of rechargeable battery technologies \cite{LLu}. Lithium-ion (Li-ion) batteries are among the most popular rechargeable batteries available in the market. They have high energy density and low self-discharge compared to other type of batteries at the same size. The main drawback is the risk of unexpected ignition from the heat generated by the battery. A slightly improvement in term of the energy density for commercial batteries is offered by lithium polymer (LiPo) batteries. LiPo batteries are common choice to power Unmanned Aerial Vehicles (UAVs). Acceleration and maneuvering over time can cause the power to drop from peak to net power. Thus, reduce the performance of the battery significantly. In an extreme case, battery failure due to overheating can cause loss of control of the UAVs leading to property damage or casualties. Thus, there is an incentive to estimate the remaining battery capacity and its performance in an UAV operation accurately. Other types of rechargeable batteries notable for their high energy density include lithium-sulfur (Li-S) and sodium-ion batteries (SIB). The former is still under development but may potentially replace lithium-ion cells not only because of its higher energy density, but also because its reduced cost from the use of sulfur \cite{George}.

The safety systems for rechargeable batteries can be improved using a smart battery management system (BMS). The BMS monitors the states and parameters of the batteries, thus provides instantaneous information about when a battery needs to be recharged and allows the BMS to prolong the battery life by preventing batteries from over-charging or over discharging.  An important component of BMS is an accurate estimation of state-of-charge (SoC) and state-of-health (SoH) \cite{KF2}. The SoC is defined as the available battery capacity, measured in Ampere-hours, and expressed as a percentage of its rated capacity. The SoC parameters can be seen as thermodynamic quantities enabling one to assess the potential energy of a battery. The SoH represents battery performance to store and deliver electrical energy relative to a new battery. The SoC and SoH cannot be measured directly and has to be estimated based on measurable variables such as current and terminal voltage. The estimation is usually done by employing mathematical or physical models to describe complex nonlinear dynamic processes arising from thermodynamics, electrode kinetics, and transport phenomena \cite{1RC}. The estimate can be used to detect battery anomaly such as overheating in failure prevention system.

There are two common approaches for SoC and SoH estimations of rechargeable batteries:
\begin{itemize}
\item \textbf{Kalman filter-based approach}. This is a natural choice since the model and the measurements in rechargeable batteries are subject to noises. Kalman filter was designed to suppress the noise affecting the battery systems. The linear Kalman filter is globally exponentially stable and optimal in the sense of minimum variance, as long as the model perfectly matches the real system, the noise is white and uncorrelated, and the covariances of the noise are known. However, battery models are nonlinear. The nonlinearities came from the measurement equation. To apply Kalman filter for SoC estimation, the system has to be linearized. The use of linearized model for the covariance error is known as the Extended Kalman filter (EKF) and has been used extensively in battery SoC estimation \cite{EKF1,EKF2,EKF3,EKF4,EKF6}. However, this method come at a price since the stability and optimality of the Kalman filter could not hold anymore. Some other variants include the Unscented Kalman filter (UKF) \cite{UKF1,UKF2,UKF3} combined with genetic algorithm \cite{Rui1,Rui2} and the particle filter \cite{Par1}. However, they are expensive in terms of computational time and are still approximations without theoretical guarantees of stability and optimality.
\item \textbf{Observer-based approach}. This approach offers global asymptotic or exponential stability. The design involves a tuning parameter that can be adjusted to pursue desired performance. Usually the analysis is based on Lyapunov stability theory. This approach include nonlinear observer \cite{OB1,OB3,OB4}, sliding-mode observer \cite{Slid}, backstepping PDE observer \cite{Moura,A4,A5}, and proportional-integral observer \cite{OB2}. These methods are relatively straightforward to implement and has higher computational efficiency.
\end{itemize}

In this paper, we present a new framework for state-of-charge estimation of rechargeable batteries using a two-stage nonlinear estimator, called the eXogenous Kalman filter (XKF), developed in \cite{Tor}. The XKF is a cascade of a nonlinear observer (NLO) and a linearized Kalman filter (LKF). The aim is to utilize each strength from both nonlinear state estimation methods, such as the global stability property of the NLO and the optimality of the LKF. The main idea is to use the exogenous state estimation obtained from the NLO to generate a linearized model for the Kalman filter.

\begin{figure}[h!]
  \centering
      \includegraphics[width=1\textwidth]{schematic}
  \caption{The schematic diagram of the XKF for SoC estimation. The XKF is a cascade of NLO and LKF.}
\label{3}
\end{figure}

In Fig. \ref{3}, we denote $\overline{\mathbf{x}}$ as the exogenous state estimation of SoC from the NLO, while $\hat{\mathbf{x}}$ denotes the SoC estimation from the LKF. The nonlinear observer gain $\mathbf{K}$ could be obtained from an analysis of the error system using Lyapunov stability theory, while the Kalman gain $\mathbf{K}_f$ is obtained after solving a Riccati equation. It has been shown that the two-stage nonlinear estimator inherits the global stability property of the NLO \cite{Tor}.

To illustrate our proposed method, we choose to work with lithium iron phosphate (LiFePO$_4$) battery. LiFePO$_4$ offers a longer life cycle and has a more stable discharge voltage than other lithium-ion batteries. It is mainly used in electric vehicles and solar energy installations \cite{Rahul}. In section 2, we present a model for LiFePO$_4$ based on Coulomb counting and 2-RC equivalent circuit model (ECM). This model is chosen because of its simple structure, thus suitable for real-time application. In section 3, we design a NLO and a LKF for our 2-RC ECM. We prove the NLO to be globally exponentially stable. Evaluation using experimental data from Dynamic Stress Test (DST) and Federal Urban Driving Schedule (FUDS) are presented in section 4. Finally, section 5 contains conclusions and future works.

\section{Lithium-ion Battery Modeling}

A lithium-ion battery uses a metal oxide cathode, a porous carbon anode, a separator, and electrolyte as conductor. The separator allows the transport of lithium-ions from one side to the other, and at the same time prevents electrons from passing through. The electrolyte is electrically conductive and contains free ions in which the lithium-ions can be transported easily \cite{John}. During discharge, the ions flow from the anode to the cathode through the electrolyte and a separator. Charging will reverse the direction and the ions flow from the cathode to the anode.

There are several types of lithium-ion batteries based on the chemistry, performance, cost, and safety. A LiFePO$_4$ battery, also called LFP battery, uses LiFePO$_4$ as a cathode material, and a graphitic carbon electrode with a metallic current collector grid as the anode. The chemical reactions in the cathode and the anode of the LiFePO$_4$ battery could be described as follow
\begin{eqnarray}
\ce{LiFePO_4 &<=>[charge][discharge]& FePO_4 + Li^+ + e^-}\nonumber\\
\ce{Li^+ + e^- + C_6 &<=>[charge][discharge]& LiC_6}\nonumber
\end{eqnarray}
The original mineral of LiFePO$_4$ is olivine which exists as a solid solution of the pure magnesium phase, forsterite, and the purely ferrous phase fayalite. So far, LiFePO$_4$ is the safest of the mainstream lithium-ion battery types. A LiFePO$_4$ battery does not require to be fully charged. Its performance even slightly improves in case of partial charge instead of a full charge. This becomes a major advantage of LiFePO$_4$ compared to lead-acid battery types. Other advantages include wide operating temperature range, excellent cycling performance, low internal resistance, and high efficiency. These are some reasons why LiFePO$_4$ has been used in electric vehicle and backup power.

A LiFePO$_4$ battery has also some drawbacks. The cell will fail if the voltage over the cell falls to less than 2.5V or increases to more than 4.2V. Furthermore, the cells do not auto-balance at the end of the charge cycle. A BMS could help to balance these cells. Furthermore, it could prevent cell under voltage by timely disconnecting the load and cell overvoltage by reducing charge current or stopping the charge process. A key component in BMS is a good model to describe the behavior of the battery under some operational condition. Lithium-ion battery modeling could be done in three different approaches:
\begin{itemize}
\item \textbf{Distributed Parameter Model}. This approach is based on thermal-electrochemical model \cite{Mou1,Hasan,Shux}. The model is derived from first principles and is written as a system of coupled partial differential equations (PDEs) and ordinary differential equations (ODEs) to describe the physics of the battery. Examples of this approach are the Doyle-Fuller-Newman model \cite{Doyle} and the single particle model (SPM) \cite{Wh}. This approach can accurately match experimental data but are complex and time consuming to solve. Advanced battery modeling using coupled PDEs-ODEs has received much attention in recent years since the estimation method for this kind of systems has been well established \cite{A1,A2,A20,A30}.
\item \textbf{Lumped Model}. This approach considers a battery as an equivalent circuit that represents the electrical characteristics of the cell \cite{Eval}. The model, which is taken from the RC ECMs, serves as a proxy model and is written as a system of ODEs. The lumped model can be solved quickly and used in real-time environments \cite{A3}. The main drawback is that this approach does not accurately capture the physics of the battery at high discharge conditions.
\item \textbf{Data-driven Model}. The last approach is based on data-driven and relies on machine learning algorithms to develop relationships between sensor inputs and output of the battery \cite{Neu1}. This approach usually requires significant amounts of hand-crafted structured training data and do not learn incrementally or interactively in real time. Next generation of machine learning technology is being researched.
\end{itemize}

In mobile platforms like unmanned aerial vehicles (UAVs), the complexity and accuracy of the model should be well balanced \cite{AThor}. Therefore, we present a lithium-ion battery model based on the ECM. An ECM uses a circuit consisting of voltage sources, resistors, and capacitors to simulate the dynamic characteristics of batteries \cite{Zhong}. The ECM offers a low complexity description of the lithium-ion battery dynamics with fewer states and parameters, thus more applicable in terms of practicality compared to the high-fidelity distributed parameter model. The battery SoC quantifies the usable energy at the present cycle and can be calculated by
\begin{eqnarray}
\text{SoC}(t) &=& \text{SoC}(t_0) - \frac{\int_{t_0}^{t} I(\tau)\mathrm{d}\tau}{Q_c}\label{main}
\end{eqnarray}
where $t_0$ denotes the initial time, $I$ is the electrical current, and $Q_c$ is the nominal capacity. Calculating the SoC from the measured discharging current and integrating it over time is known as Coulomb counting. This method, however, suffers from long-term drift and lack of a reference point. Therefore, the SoC should be calibrated on regular basis. Measured terminal voltages can be used for calibration. The terminal voltage can be obtained from the open circuit voltage (OCV), which underlying physical phenomenon of lithium-ion intercalation/deintercalation process and is expressed as a nonlinear function of the SoC. In practice, we create an OCV-SoC curve from experiments and use it either as an analytical expression or as a lookup table for the model. Expressing the OCV-SoC curve analytically using polynomial function is prone to incorrect estimation, since this approach is very sensitive with respect to the polynomial coefficient. Therefore, we chose to use the OCV-SoC as a lookup table. To calculate its derivative for the LKF, we approximate it using finite difference method from its interpolation. An example of an OCV-SoC curve for LiFePO$_4$ batteries could be seen in Fig. \ref{1}. The OCV-SoC curve is a function of the ambient temperature. Note that a small inaccuracy in the OCV-SoC curve may cause a large deviation in the SoC estimation, since there is a wide flat region on the OCV-SoC curve \cite{Weng}. Thus, improving the OCV models can contribute significantly to increase the accuracy of SoC estimation.
\begin{figure}[h!]
  \centering
      \includegraphics[width=1\textwidth]{SOCTy}
  \caption{Typical OCV-SoC lookup tables for LiFePO$_4$ battery. The ambient temperature alters the OCV-SoC lookup tables.}
\label{1}
\end{figure}

Using the OCV-SoC lookup table and the ECM (Fig. \ref{2}), the terminal voltage $V$ can be calculated using the following formula
\begin{eqnarray}
V(t) &=& V_{oc}(\text{SoC})-I(t)R-C(t)\label{main2}
\end{eqnarray}
where $R$ is the internal resistance or the ohmic resistance. This ohmic resistance is used to represent the electrical resistance of battery components with the accumulation and dissipation of charge in the electrical double-layer. It could be easily observed from Fig. \ref{1} that $\frac{\mathrm{d} V_{oc}}{\mathrm{d}\text{SoC}}>0$. Here, $C$ is a correction factor due to model inaccuracy and environmental conditions, e.g., ambient temperature variations. Some authors use a resistor model and consider the value of $C$ as a function of ambient temperature \cite{MP1}, which could be determined using least-square fitting from experimental data. Other authors consider $C$ as a 1-RC circuit model \cite{Xhu} and 2-RC circuit model \cite{XChen}. It could be also determined using a combination of Thevenin-based ECM \cite{tuj} with the hysteresis voltage dynamics \cite{EKF5}, which offers a grasp of dynamic current-voltage characteristics and compensates the static current-voltage property.
\begin{figure}[h!]
  \centering
      \includegraphics[width=1\textwidth]{RC}
  \caption{Schematic diagram of the 2-RC circuit for LiFePO$_4$ battery.}
  \label{2}
\end{figure}

To balance between the computational effort and accuracy, in this paper we use the 2-RC circuit model, i.e., $C(t)=V_1(t)+V_2(t)$, where
\begin{eqnarray}
\frac{\mathrm{d}V_1}{\mathrm{d}t}(t) &=& -\frac{V_1(t)}{R_1C_1}+\frac{I(t)}{C_1}\label{main3}\\
\frac{\mathrm{d}V_2}{\mathrm{d}t}(t) &=& -\frac{V_2(t)}{R_2C_2}+\frac{I(t)}{C_2}\label{main4}
\end{eqnarray}
where $R_1,R_2$ and $C_1,C_2$ are diffusion resistance and diffusion capacitance for the RC network, respectively. These parameters could also be interpreted as the mass transport effects and dynamic voltage performance. These parameters together with the ohmic resistance could be determined from the exponential-function fitting method or a simple least-square algorithm \cite{Waag}.

\section{SoC Estimation using eXogenous Kalman Filter}

Given the nonlinear battery model and measurement equations \eqref{main}-\eqref{main4}, the SoC estimation can be formulated as a nonlinear state estimation problem.

\subsection{State space model}

Differentiating \eqref{main} with respect to $t$, we have
\begin{eqnarray}
\frac{\mathrm{d}\text{SoC}}{\mathrm{d}t}(t) &=& -\frac{I(t)}{Q_c}\label{main1}
\end{eqnarray}
From \eqref{main2}-\eqref{main1}, the state-space model for the lithium-ion batteries is given by
\begin{eqnarray}
\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}V_1(t)\\V_2(t)\\\text{SoC}(t)\end{pmatrix}&=&
\begin{pmatrix}-\frac{1}{R_1C_1}&0&0\\0&-\frac{1}{R_2C_2}&0\\0&0&0\end{pmatrix}\begin{pmatrix}V_1(t)\\V_2(t)\\\text{SoC}(t)\end{pmatrix}+\begin{pmatrix}\frac{1}{C_1}\\\frac{1}{C_2}\\-\frac{1}{Q_c}\end{pmatrix}I(t)\label{ss1}\\
V(t) &=& V_{oc}(\text{SoC})-V_1(t)-V_2(t)-RI(t)\label{ss2}
\end{eqnarray}
To simplify the presentation, we write \eqref{ss1}-\eqref{ss2} as
\begin{eqnarray}
\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}(t) &=& \mathbf{A}\mathbf{x}(t)+\mathbf{B}I(t)\label{mi1}\\
V(t) &=& h(\mathbf{x})-RI(t)\label{mi2}
\end{eqnarray}
where
\begin{eqnarray}
\mathbf{x} &=& \begin{pmatrix}V_1(t)\\V_2(t)\\\text{SoC}(t)\end{pmatrix}\\
\mathbf{A} &=& \begin{pmatrix}-\frac{1}{R_1C_1}&0&0\\0&-\frac{1}{R_2C_2}&0\\0&0&0\end{pmatrix}\\
\mathbf{B} &=& \begin{pmatrix}\frac{1}{C_1}\\\frac{1}{C_2}\\-\frac{1}{Q_c}\end{pmatrix}\\
h(\mathbf{x}) &=& V_{oc}(\text{SoC})-V_1(t)-V_2(t)
\end{eqnarray}
To incorporate the model and measurement uncertainties and inaccuracies, noises are added into \eqref{mi1}-\eqref{mi2}, thus the complete model becomes
\begin{eqnarray}
\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}(t) &=& \mathbf{A}\mathbf{x}(t)+\mathbf{B}I(t)+w(t)\\
V(t) &=& h(\mathbf{x})-RI(t)+v(t)\\
w &\sim& (0,\mathbf{Q}_n)\\
v &\sim& (0,\mathbf{R}_n)
\end{eqnarray}
where $w$ denotes the process noise, $\mathbf{Q}_n$ is the process noise covariance, $v$ is the measurement noise, and $\mathbf{R}_n$ is the measurement noise covariance. To use the optimal filter, the covariances of the noise $\mathbf{Q}_n$ and $\mathbf{R}_n$ have to be known. These values can be obtained using Bayesian, maximum likelihood, covariance matching, and correlation techniques \cite{Row}. Furthermore, the linear Kalman filter is optimal under the assumption that $w$ and $v$ are white noises and uncorrelated. In this section, we design the SoC estimation using the XKF, which is a cascade of a NLO and LKF.

\subsection{Nonlinear Observer (NLO)}

The NLO is designed as follow
\begin{eqnarray}
\frac{\mathrm{d}\overline{\mathbf{x}}}{\mathrm{d}t}(t)(t) &=& \mathbf{A}\overline{\mathbf{x}}(t)+\mathbf{B}I(t)+\mathbf{K}\left(V(t)-\overline{V}(t)\right)\label{no1}\\
\overline{V}(t) &=& h(\overline{\mathbf{x}})-RI(t)\label{no2}
\end{eqnarray}
where $\overline{\mathbf{x}}$ denotes the exogenous state estimation from the NLO and $\mathbf{K}=\begin{pmatrix}k_1&k_2&k_3\end{pmatrix}^{\intercal}\in\mathbb{R}^3$ is the observer gain to be determined later. Let us define the difference between the actual state and the measured state as $\breve{\mathbf{x}}(t) = \mathbf{x}(t) - \overline{\mathbf{x}}(t)$ and the difference between the actual terminal voltage and the measured terminal voltage as $\breve{V}(t) = V(t)-\overline{V}(t)$. The error equation is obtained as follow
\begin{eqnarray}
\frac{\mathrm{d}\breve{\mathbf{x}}}{\mathrm{d}t}(t) &=& \mathbf{A}\breve{\mathbf{x}}(t)-\mathbf{K}\breve{V}(t)\label{err}\\
\breve{V}(t) &=& h(\mathbf{x}) - h(\overline{\mathbf{x}})\label{meas}
\end{eqnarray}
Utilizing the smoothness property of the OCV-SoC lookup table, we can linearize the nonlinear measurement equation \eqref{meas} using the mean value theorem. Note that the continuity of the OCV-SoC lookup table can be obtained by applying interpolation from the experimented data. First, we write
\begin{eqnarray}
\breve{V}(t) &=& -\breve{V}_1(t)-\breve{V}_2(t)+V_{oc}(\text{SoC})-V_{oc}(\overline{\text{SoC}})
\end{eqnarray}
From the mean value theorem, there exist $\text{SoC}\leq\xi\leq\overline{\text{SoC}}$, such that
\begin{eqnarray}
V_{oc}(\text{SoC})-V_{oc}(\overline{\text{SoC}}) = \frac{\mathrm{d}V_{oc}}{\mathrm{d}\text{SoC}}(\xi)\breve{\text{SoC}}(t)
\end{eqnarray}
where $\breve{\text{SoC}}(t)=\text{SoC}(t)-\overline{\text{SoC}}(t)$. Let us denote
\begin{eqnarray}
\mathbf{T} &=& \begin{pmatrix}-1&-1&\frac{\mathrm{d}V_{oc}}{\mathrm{d}\text{SoC}}(\xi)\end{pmatrix}\label{TT}
\end{eqnarray}
Thus, we can write
\begin{eqnarray}
\breve{V}(t) = \mathbf{T}\breve{\mathbf{x}}(t)
\end{eqnarray}
and the error equation becomes
\begin{eqnarray}
\frac{\mathrm{d}\breve{\mathbf{x}}}{\mathrm{d}t}(t) &=& \left(\mathbf{A}-\mathbf{K}\mathbf{T}\right)\breve{\mathbf{x}}(t)
\end{eqnarray}
Note that $\frac{\mathrm{d}V_{oc}}{\mathrm{d}\text{SoC}}(\xi)>0$ and $\mathbf{T}$ is time-varying and bounded. In this particular case, let $k_1=k_2=0$. Then, $\breve{V}_1(t)=\breve{V}_1(0)e^{-\frac{1}{R_1C_1}}$ and $\breve{V}_2(t)=\breve{V}_2(0)e^{-\frac{1}{R_2C_2}}$, i.e., $\breve{V}_1(t)$ and $\breve{V}_2(t)$ are exponentially stable. Correspondingly, we have $\dot{\breve{\text{SoC}}}(t)=-k_3\frac{\mathrm{d}V_{oc}}{\mathrm{d}\text{SoC}}(\xi)\breve{\text{SoC}}(t)$. Thus, if $k_3>0$ then the nonlinear observer is globally exponentially stable, i.e., the error will decay to zero and the estimate $\overline{\mathbf{x}}$ will converge to the actual value exponentially.

\subsection{Linearized Kalman Filter (LKF)}

Now we have $\overline{\mathbf{x}}$ as an estimate of $\mathbf{x}$, which is a bounded signal given by the NLO \eqref{no1}-\eqref{no2}. We use this signal as a linearization point for the linear Kalman filter. A first-order Taylor series expansion about the trajectory $\overline{\mathbf{x}}$ gives the linearized model
\begin{eqnarray}
\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}(t) &=& \mathbf{A}\mathbf{x}(t)+\mathbf{B}I(t)+w(t)\\
V(t) &=& h(\overline{\mathbf{x}})+\mathbf{H}(\overline{\mathbf{x}})\breve{\mathbf{x}}+r(\mathbf{x},\overline{\mathbf{x}})-RI(t)+v(t)
\end{eqnarray}
where
\begin{eqnarray}
\mathbf{H}(\overline{\mathbf{x}}) &=& \frac{\partial h}{\partial\mathbf{x}}(\overline{\mathbf{x}}) = \begin{pmatrix}-1&-1&\frac{\mathrm{d} V_{oc}}{\mathrm{d}\text{SoC}}(\overline{\text{SoC}})\end{pmatrix}
\end{eqnarray}
Remark that, since $\overline{\mathbf{x}}$ is bounded and converges to $\mathbf{x}$, we can neglect the higher-order term $r(\mathbf{x},\overline{\mathbf{x}})$ in the design since it has no consequences for stability \cite{Tor}. Thus, we can design the second stage estimator $\hat{\mathbf{x}}$ using the LKF as follows
\begin{eqnarray}
\frac{\mathrm{d}\hat{\mathbf{x}}}{\mathrm{d}t}(t) &=& \mathbf{A}\hat{\mathbf{x}}(t)+\mathbf{B}I(t)+\mathbf{K}_f(t)\left(V(t)-h(\overline{\mathbf{x}})-H(\overline{\mathbf{x}})(\hat{\mathbf{x}}-\overline{\mathbf{x}})+RI(t)\right)\label{xkf}
\end{eqnarray}
The time-varying gain $\mathbf{K}_f(t)$ is given by
\begin{eqnarray}
\mathbf{K}_f(t) = \mathbf{P}(t)\mathbf{H}^{\intercal}(\overline{\mathbf{x}})\mathbf{R}_n^{-1}
\end{eqnarray}
where the error covariance matrix $\mathbf{P}(t)$ is the solution to the Riccati equation
\begin{eqnarray}
\frac{\mathrm{d}\mathbf{P}}{\mathrm{d}t}(t) &=& \mathbf{A}\mathbf{P}(t)+\mathbf{P}(t)\mathbf{A}^{\intercal}+\mathbf{Q}_n-\mathbf{K}_f(t)\mathbf{R}_n\mathbf{K}_f^{\intercal}(t)
\end{eqnarray}
with $\mathbf{P}(0)$ symmetric and positive definite. Note that, unlike EKF, the XKF uses linear time-varying measurement model that is independent of the estimate $\hat{\mathbf{x}}$ in \eqref{xkf}. To implement the algorithm with a discrete-time Kalman filter, the model has to be discretized, for example using the Euler discretization method.

\section{Evaluation using Experimental Data}

The LiFePO$_4$ battery is tested in two dynamic loading condition tests: the dynamic stress test (DST) and the federal urban driving schedule (FUDS). The DST was used to identify the model parameters, while the FUDS was used to validate the performance of the SoC estimation. The battery specification is given in Table 1. A complete description regarding the tests is given in \cite{MP1}.
\begin{table}[h!]
\begin{center}
  \begin{tabular}{ l | c | r | r}
    \hline\hline
    Type & Nominal & Nominal & Upper and lower \\
    &voltage & capacity & cut-off voltage\\ \hline
    LiFePO$_4$ & 3.3V & 2.23Ah & 3.6V and 2.0V\\
    \hline
  \end{tabular}
  \caption{Specification of the LiFePO$_4$ battery used in the experiment.}
\end{center}
\end{table}

\subsection{Dynamic Stress Test}

To evaluate the validity and to identify the parameters of the battery model, the DST is performed at 20$^\circ$C. This test is designed by US Advanced Battery Consortium (USABC) to simulate a variable-power discharge regime that represents the expected demands of an electric vehicle (EV) battery. The voltage and current are measured and recorded from fully charged to empty with a sampling period of 1.0s based on the battery test bench. The accumulative charge was run continuously from 100\% SoC at 3.6V to empty at 2.0V over several cycles in a discharge process. The measured current $I$ and voltage $V$ are given in Fig. \ref{4}. This measurements are used to calibrate the estimation from the Coulomb counting method.
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{MeaCurrent}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{MeaVoltage}
\end{subfigure}
\caption{Measured current and terminal voltage under the dynamic stress test (DST).}
\label{4}
\end{figure}

The root mean square (RMS) error is used to evaluate the validity of the model and the parameters and is given by
\begin{eqnarray}
\text{RMS} = \sqrt{\frac{1}{n}\Sigma_{i=1}^n\left(V_{\text{measured}}-V_{\text{estimated}}\right)^2}
\end{eqnarray}

The model parameters and the RMS error are given in Table 2. Notice that the diffusion resistance and the diffusion capacitance for $V_1$ and $V_2$ are the same. The value of $C_1$ and $C_2$ are dominant, which means the 2-RC model gave a minor effect to SoC estimation. The accuracy of the estimation is mostly determined by the OCV-SoC look up table, which its accurate value is hard to get. This becomes a main drawback when using the ECM despite of its simplicity.
\begin{table}[h!]
\begin{center}
  \begin{tabular}{ l | c | c | c | c | c}
    \hline\hline
    $R$ & $R_1$ & $C_1$ & $R_2$ & $C_2$ & \text{RMS}\\\hline
    0.18 & 0.035 & 1e6 & 0.035 & 1e6 & 8.4157e-05\\
    \hline
  \end{tabular}
  \caption{Model parameters and the RMS.}
\end{center}
\end{table}

Fig. \ref{5} shows the fitting between the measured and estimated terminal voltage $V$ and the error. A slightly large deviation can be observed when the voltage is close to 2.0V, or when the SoC approaches zero. This is due to inaccuracies when measuring the open circuit voltage. The OCV-SoC lookup table and the estimated SoC are given in Fig. \ref{6}. The estimated SoC is used as a reference for the XKF.
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{Voltage}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{MeasurementError}
\end{subfigure}
\caption{Measured and estimated voltage response under the DST and the model error.}
\label{5}
\end{figure}

\begin{figure}[h!]
\centering
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{OCV}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{SoC}
\end{subfigure}
\caption{The open circuit voltage and SoC under the DST.}
\label{6}
\end{figure}

To test the proposed method, we run the estimations with initial SoC at 60$\%$ and compare the results with some existing filter, e.g, EKF and UKF. The standard deviation for the process noise $\mathbf{Q}_n$ is 0.01, while for the measurement noise $\mathbf{R}_n$ is 0.04. These values are obtained using a simple covariance matching technique. The initial error covariance matrix $\mathbf{P}(0)$ is defined based on the initialization error. In this case, if the initial state is not very close, the value of $\mathbf{P}(0)$ should be large, whereas if the initialization is very good a smaller $\mathbf{P}(0)$ value can be used. The results are given in Fig. \ref{7}. It was observed that the NLO performs better in terms of the convergence rate than the EKF and UKF. The XKF further improves the estimation. Another experiment is done by starting the estimation when the SoC is at 78\%. Using XKF at different initial estimations, the values converge to the actual SoC, as can be seen from Fig. \ref{71}.
\begin{figure}[h!]
  \centering
      \includegraphics[width=1\textwidth]{SoCDST20}
  \caption{SoC estimation from different methods under DST with true initial condition 100\%.}
\label{7}
\end{figure}

\begin{figure}[h!]
  \centering
      \includegraphics[width=1\textwidth]{SOCDST70}
  \caption{SoC estimation using XKF under DST with true initial condition 78\%.}
\label{71}
\end{figure}

The convergence time can be increased using higher value of the gain $k_3$, as can be seen from Fig. \ref{hg}. With $k_3=2$ the estimate converges to the actual value in one minutes, while with $k_3=0.5$, the estimate converges after almost one hour. However, keep in mind that higher observer tends to overshoot as can be seen between $t=6500s$ to the end. Indeed, selecting an appropriate gain in crucial in SoC estimation. There is a trade-off between convergence rate and accuracy.
\begin{figure}[h!]
  \centering
      \includegraphics[width=1\textwidth]{HighGain}
  \caption{Comparison of different observer gain $k_3$.}
\label{hg}
\end{figure}

\subsection{Federal Urban Driving Schedule}

A more realistic dynamic current test is given by the FUDS test. FUDS is based on a time-velocity profile from an automobile industry standard vehicle to test the dynamic electric vehicle performance. The measured current and voltage can be seen in Fig. \ref{8}. The model fitting and the error can be seen in Fig. \ref{10}, while the OCV-SoC lookup table and the SoC could be seen in Fig. \ref{11}, respectively. The current causes variation of the SoC from fully charged at 3.6V to empty at 2.0V. We run the simulations from three different initial guesses using XKF. The results are given in Fig. \ref{9}. It can be observed that the estimation converges to the SoC for any initial guess. The convergence rate can be increased using higher observer gain, but it  may cause significant overshoot and lead to oscillation and instability.
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{CurrentFUDS}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{VoltageFUDS}
\end{subfigure}
\caption{Measured current and terminal voltage response under the federal urban driving schedule (FUDS).}
\label{8}
\end{figure}

\begin{figure}[h!]
\centering
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{Voltage40}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{MeasurementError40}
\end{subfigure}
\caption{Measured and estimated voltage response under the FUDS and the model error.}
\label{10}
\end{figure}

\begin{figure}[h!]
\centering
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{OCV40}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{SoC40}
\end{subfigure}
\caption{The open circuit voltage and SoC under the FUDS.}
\label{11}
\end{figure}

\begin{figure}[h!]
  \centering
      \includegraphics[width=1\textwidth]{SoCFUDS}
  \caption{The SoC estimation from XKF under the FUDS at 20$^\circ$C.}
\label{9}
\end{figure}

\subsection{Effect on Different Ambient Temperatures}

In the last experiment, we estimate the SoC for different ambient temperatures. As we have discussed earlier, the OCV-SoC curve for LiFePO$_4$ has a wide flat region. Different ambient temperature will alter this curve significantly, especially at lower SoC. This can be seen from Fig. \ref{14}. The blue line shows the SoC at 20$^\circ$C, while the green line shows the SoC at 40$^\circ$C. We observe that higher ambient temperature will reduce the SoC. Furthermore, if the SoC estimation at 40$^\circ$C is calculated using the OCV-SoC curve at 20$^\circ$C (red line), then there is a deviation especially after the SoC is less than 20\%. While if it uses the OCV-SoC curve at 40$^\circ$C (black line), the estimate converges to the actual value. Thus, the OCV-SoC curve needs to be adjusted accordingly. Fail to update the OCV-SoC curve can cause incorrect SoC estimation. Thus, there is an incentive to model the temperature variation for the OCV-SoC curve, e.g, using a thermal-electrochemical model, which will be investigated in our further works.
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{TambOCV}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
   \includegraphics[width=1\linewidth]{TambSOC}
\end{subfigure}
\caption{The open circuit voltage and SoC estimation using XKF under the FUDS at 40$^\circ$C.}
\label{14}
\end{figure}

\section{Conclusions and Future Works}

In this paper, we present SoC estimation of a LiFePO$_4$ battery using the eXogenous Kalman filter from the Coulomb counting and equivalent circuit model. The estimation algorithm consists of a cascade of NLO and LKF. Simulations against experimental data from the dynamics stress test and federal urban driving schedule show the algorithm is able to estimate the SoC accurately. Furthermore, comparisons with existing filter show the estimation using XKF converges faster, thanks to the exponential stability from the NLO. Future work includes the use of thermal-electrochemical model to improve the lithium battery estimation.

\section*{Acknowledgements}

This work was supported by the Free the Drones (FreeD) project at Center for Unmanned Aerial Vehicles, University of Southern Denmark, project number 5156-00008B, and the FRINATEK project "Multi-stage global sensor fusion for navigation" and the Centre of Autonomous Marine Operations and Systems (NTNU AMOS) at the Norwegian University of Science and Technology (NTNU), project numbers 250725 and 223254, respectively.