Liang Li

Stéphane Lanteri
email: stephane.lanteri@inria.fr

N Asger Mortensen

Martijn Wubs
email: mwubs@fotonik.dtu.dk

A hybridizable discontinuous Galerkin method for solving nonlocal optical response models

Keywords: Maxwell's equations, nonlocal hydrodynamic Drude model, general nonlocal optical response theory, hybridizable discontinuous Galerkin method

We propose Hybridizable Discontinuous Galerkin (HDG) methods for solving the frequencydomain Maxwell's equations coupled to the Nonlocal Hydrodynamic Drude (NHD) and Generalized Nonlocal Optical Response (GNOR) models, which are employed to describe the optical properties of nano-plasmonic scatterers and waveguides. Brief derivations for both the NHD model and the GNOR model are presented. The formulations of the HDG method are given, in which we introduce two hybrid variables living only on the skeleton of the mesh. The local field solutions are expressed in terms of the hybrid variables in each element. Two conservativity conditions are globally enforced to make the problem solvable and to guarantee the continuity of the tangential component of the electric field and the normal component of the current density. Numerical results show that the proposed HDG methods converge at optimal rate. We benchmark our implementation and demonstrate that the HDG method has the potential to solve complex nanophotonic problems.

Introduction

Nanophotonics is the active research field field concerned with the study of interactions between nanometer scale structures/media and light, including near-infrared, visible, and ultraviolet light. It bridges the micro and the macro worlds, and there are many connections between theoretical studies and feasible engineering. The many fascinating (potential) applications include invisibility cloaking, nano antennas, metamaterials, novel biological detection and treatment technologies, as well as new storage media [START_REF] Maier | Plasmonics: Fundamentals and Applications[END_REF].

All of the above applications of nanophotonics require elaborate control of the propagation of light waves. In order to do so, appropriate mathematical models are needed to predict the behavior of light-matter interactions. Metals are interesting for nanophotonics because they can both enhance and confine optical fields, making plasmonics of interest to emerging quantum technologies [START_REF] Tame | Quantum Plasmonics[END_REF][START_REF] Bozhevolnyi | Plasmonics for emerging quantum technologies[END_REF][START_REF] Fitzgerald | Quantum Plasmonics[END_REF]. This is enabled by the existence of Surface Plasmons (SPs). SPs are coherent oscillations that exist as evanescent waves at both sides of the interface between any two materials where the real part of the dielectric function changes sign across the interface. The typical example is a metal-dielectric interface, such as a metal sheet in air [START_REF] Sarid | Modern Introduction to Surface Plasmons: Theory, Mathematica Modeling, and Applications[END_REF]. Maxwell's equations can be employed to model the macroscale electromagnetic waves and armed with classical electrodynamics there are numerous approaches ranging from classical electrodynamics to ab initio treatments [START_REF] Gallinet | Numerical methods for nanophotonics: standard problems and future challenges Laser[END_REF][START_REF] Varas | Quantum plasmonics: from jellium models to ab initio calculations[END_REF]. Ab initio techniques can be used to simulate the microscopic dynamics on the atomic scale, but with ab initio methods one can only deal with systems with up to about ten thousand atoms [START_REF] Varas | Quantum plasmonics: from jellium models to ab initio calculations[END_REF], thus calling for semiclassical treatments [START_REF] Mortensen | A generalized non-local optical response theory for plasmonic nanostructures[END_REF][START_REF] Toscano | Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics[END_REF] or more effective inclusions of quantum phenomena into classical electrodynamics [START_REF] Luo | Surface plasmons and nonlocality: a simple model[END_REF][START_REF] Yan | Projected dipole model for quantum plasmonics[END_REF][START_REF] Christensen | Quantum corrections in nanoplasmonics: shape, scale, and material[END_REF][START_REF] Zhu | Quantum mechanical effects in plasmonic structures with subnanometre gaps[END_REF].

If one models the interaction of light with metallic nanostructures classically or semiclassically, then this calls for appropriate modelling of the material response as described for example by the Drude model [START_REF] Drude | Zur Elektronentheorie der Metalle[END_REF][START_REF] Dressel | Verifying the Drude response[END_REF], the Nonlocal Hydrodynamic Drude (NHD) model [START_REF] Bloch | Bremsvermögen von Atomen mit mehreren Elektronen[END_REF][START_REF] Raza | Nonlocal optical response in metallic nanostructures[END_REF], or the Generalized Nonlocal Optical Response (GNOR) theory [START_REF] Mortensen | A generalized non-local optical response theory for plasmonic nanostructures[END_REF], all in combination with and coupled to Maxwell's equations. Except for some highly symmetric geometries, analytical solutions to the resulting systems of differential equations are not available. Thus, numerical treatment of these systems of PDEs is an important aspect of nanophotonics research. Numerical experiments help to find promising systems and geometries before real fabrication, to obtain optimized parameters, to visualize field distributions, to investigate the dominant contribution to a phenomenon, to explain experimental observations, and so on [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF].

Several numerical methods exist for computing the solution of Maxwell's equations [START_REF] Gallinet | Numerical methods for nanophotonics: standard problems and future challenges Laser[END_REF]. For time-dependent problems, the Finite-Difference Time-Domain (FDTD) algorithm is the most popular method [START_REF] Taflove | Computational Electrodynamics -The Finite-Difference Time-Domain Method[END_REF] among physicists and engineers. More recently, the Discontinuous Galerkin Time-Domain (DGTD) method has drawn a lot of attention because of several appealing features, for example, easy adaptation to complex geometries and material composition, high-order accuracy, and natural parallelism [START_REF] Hesthaven | Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell's equations[END_REF]. For time-harmonic problems, the Finite Element Method (FEM) is most widely used for the solution of Maxwell's equations. In very recent years, the Hybridizable Discontinuous Galerkin (HDG) method appears as a promising numerical method for time-harmonic problems because it inherits nearly all the advantages of the DG methods while leading to a computational complexity similar to FEM [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF][START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF][START_REF] Li | A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell's equation[END_REF][START_REF] He | Optimized Schwarz algorithms for solving time-harmonic Maxwell's equations discretized by a hybridizable discontinuous Galerkin method[END_REF].

Currently, FDTD (for time-dependent problems) and FEM (for time-harmonic problems) methods are still the methods most commonly adopted for the simulation of light-matter

Incident plane wave

Figure 1: Sketch of the incident electromagnetic wave illuminating the scatterer Ω S that has a subwavelength size and is surrounded by free space. Ω S is usually filled with metals, such as gold, silver or sodium. An artificial absorbing boundary ∂Ω is introduced to make a computational domain Ω.

interactions. Most often, commercial simulation software (such as Lumerical FDTD1 and Comsol Multiphysics2 ) is used for that purpose. However, these methods and computer codes do not always offer the required capabilities for addressing accurately and efficiently the complexity of the physical phenomena underlying nanometer scale light-matter interactions. In the academic community, also the DGTD method has recently been considered in this context [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF][START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF][START_REF] Huang | Theoretical and numerical analysis of a non-local dispersion model for light interaction with metallic nanostructures[END_REF]. In Ref. [START_REF] Schmitt | A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects[END_REF], some numerical results are presented for the NHD model using the DGTD method. In the present paper we are employing the HDG method to solve the frequency-domain NHD and GNOR models. The development of accurate and efficient numerical methods for computational nanophotonics is expected to be a long-lasting demand, both because new models are regularly proposed that require innovative numerical methods, and because there is demand for more accurate and faster simulation methods for existing models.

This paper introduces a HDG method for the solution the NHD and GNOR models. The rest of the paper is organized as follows. In section 2, we briefly introduce mathematical aspects both of the NHD model and of the GNOR model. HDG formulations are given in section 3. Numerical results are presented in section 4 to show the effectiveness of high-order HDG methods for solving problems in nanophotonics. We draw conclusions in section 5.

Physics problem: nonlocal optical response by nanoparticles

The problem considered is shown in Figure 1 where the nanometer-size metal Ω S is illuminated by an incident plane wave of light. The infinite scattering domain is truncated as a finite computational domain Ω by employing an artificial absorbing boundary condition, which is designed to absorb outgoing waves.

Nonlocal hydrodynamic Drude model

There are a number of theories for the modeling of the light-matter interactions which are used under different settings. In this subsection, we briefly introduce the NHD model. The incoming light is described as a propagating electromagnetic wave that satisfies Maxwell's equations. Without external charge and current, Maxwell's equations of macroscopic electromagnetism for non-magnetic materials can be written as

∇ × H = ε 0 ε loc ∂ t E + J, ∇ × E = -µ 0 ∂ t H, (1)
where H and E are respectively the magnetic and electric fields, ε 0 is the permittivity constant, µ 0 is the permeability constant, ε loc = ε ∞ + ε inter is introduced to account for the local response, and J is the nonlocal hydrodynamic polarization current density which is due to the nonlocal material on the plasmonic scatterers [START_REF] Hiremath | Numerical solution of nonlocal hydrodynamic Drude model for arbitrary shaped nano-plasmonic structures using Nédélec finite elements[END_REF]. In this paper, we will for simplicity set ε inter = 0 and ε ∞ = 1, thereby focusing solely on the free-electron response to light. Equations (1) need to be completed to solve electromagnetic fields E and H because of the unknown polarization current density J. The models that we will consider in this paper differ only in the assumed dynamics of the polarization current density, which we will now discuss in more detail. The polarization current density J due to the motion of the free-electron gas can be written as

J = -env, (2)
where e is the charge of the electron, n is the density of the electron gas (a scalar field), and v is its hydrodynamic velocity (a vector field). Within the hydrodynamic model, the dynamics of the velocity field is given by [START_REF] Schmitt | A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects[END_REF][START_REF] Raza | Nonlocal study of ultimate plasmon hybridization[END_REF]

m e (∂ t + v • ∇)v = -e(E + v × B) -m e γv -∇ δg[n] δn , (3)
where m e is the mass of an electron, -e(E + v × B) is the Lorentz force with B being the magnetic flux density, γ is a damping constant, g[n] is an energy functional of the fluid, and the term ∇ δg [n] δn denotes the quantum pressure. Complementary to Eq. 3, the dynamics of the free-electron density is given by

∂ t n + ∇ • (nv) = 0, (4)
which is the well-known continuity relation that relates the velocity v and the density n.

The hydrodynamic dynamics described by Eq. ( 3) is obviously nonlinear in v, but in the following we only consider the linear response of the electron gas on external fields. One can write a perturbation expansion v v 0 + v 1 and similarly for the electric and magnetic fields and for the density. Since in the absence of an external field v = v 0 = 0, both the nonlinear term v • ∇v and the magnetic induction field B disappear due to the linearization [START_REF] Toscano | Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics[END_REF]. If we furthermore assume the energy functional to be of the Thomas-Fermi form, then we obtain for the linearized quantum pressure

-∇ δg[n] δn = -m e β 2 1 n 0 ∇n, (5)
where

β 2 = 3 5 v 2
F with v F being the Fermi velocity. The zero-order (i.e. equilibrium) density n 0 is constant within the plasmonic medium [START_REF] Toscano | Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics[END_REF]. Here in Eq. ( 5) and below, we write n for the linearized density n 1 and similarly we will from now on simply write v for the linearized velocity v 1 . As a result, we obtain the linearized hydrodynamic equation [START_REF] Raza | Nonlocal optical response in metallic nanostructures[END_REF][START_REF] Schmitt | A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects[END_REF]

m e ∂ t v = -eE -m e γv -m e β 2 1 n 0 ∇n, (6)
as well as the linearized continuity relation

∂ t n = -n 0 ∇ • v. (7)
Inserting Eqs. (2) (linearized as J = -en 0 v) and ( 7) into [START_REF] Gallinet | Numerical methods for nanophotonics: standard problems and future challenges Laser[END_REF], and taking the time-derivative ∂ t , we obtain

∂ tt J + γ∂ t J -β 2 ∇(∇ • J) -ω 2 p ε 0 ∂ t E = 0, (8)
where ω p is the plasma frequency with ω 2 p = n 0 e 2 /(m e ε 0 ). By Fourier transformation we replace ∂ t with -iω, where i is the imaginary unit and ω is the angular frequency, and obtain the frequency-domain relation between polarization current density and the electric field within the hydrodynamic model as

ω(ω + iγ)J + β 2 ∇(∇ • J) = iωω 2 p ε 0 E. (9)
This equation describes electron-field interaction within the plasmonic nanostructure Ω S . We will neglect spill-out of electrons outside the classical geometric surface of the structure, which for our purposes is a good assumption for noble metals such as silver and gold [START_REF] Toscano | Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics[END_REF]. Mathematically, this is arranged by imposing a hard-wall condition on the boundary ∂Ω S , namely n • J = 0 on ∂Ω S [START_REF] Jewsbury | Electrodynamic boundary conditions at metal interfaces[END_REF][START_REF] Yan | Green's function surface-integral method for nonlocal response of plasmonic nanowires in arbitrary dielectric environments[END_REF].

General nonlocal optical response model

We also briefly present the mathematical derivation of the central equations of the GNOR model, based on Ref. [START_REF] Mortensen | A generalized non-local optical response theory for plasmonic nanostructures[END_REF]. In the GNOR model, also diffusion of the electron gas is taken into consideration. Let the density n(r, t) = n 0 + n 1 (r, t), where the last term is the induced density variation caused by a non-vanishing electric field E, which we assume sufficiently small that n 1 n 0 holds. Instead of [START_REF] Varas | Quantum plasmonics: from jellium models to ab initio calculations[END_REF], we now consider the linearized convection-diffusion equation [START_REF] Mortensen | A generalized non-local optical response theory for plasmonic nanostructures[END_REF] ∂

t (-e 1 )n 1 = D∇ 2 (-e)n 1 -∇ • [(-e)n 0 v] = -∇ • J, ( 10
)
where D is the diffusion constant for the charge-carrier diffusion. Then the current density is given by Fick's law

J = (-e)n 0 v -D∇(-e)n 1 . (11)
Multiplying ( 6) by the charge of the electron -e, the equilibrium density n 0 and taking the time-derivative we have

m e (∂ t + γ)∂ t [(-e)n 0 v] = n 0 e 2 ∂ t E -m e β 2 ∇[∂ t (-e)n 1 ]. (12)
Dividing [START_REF] Christensen | Quantum corrections in nanoplasmonics: shape, scale, and material[END_REF] by m e and combining with Fick's law [START_REF] Yan | Projected dipole model for quantum plasmonics[END_REF] results in

(∂ t + γ){∂ t J + D∇[∂ t (-e)n 1 ]} = n 0 e 2 m e ∂ t E -β 2 ∇[∂ t (-e)n 1 ]. ( 13
)
From the convection-diffusion equation ( 10), we have

(∂ t + γ)[∂ t J + D∇(∇ • J)] = n 0 e 2 m e ∂ t E -β 2 ∇(∇ • J). (14)
Like what we did for [START_REF] Mortensen | A generalized non-local optical response theory for plasmonic nanostructures[END_REF], transforming [START_REF] Drude | Zur Elektronentheorie der Metalle[END_REF] to the frequency domain gives

ω(ω + iγ)J + [β 2 + D(γ -iω)]∇(∇ • J) = iωω 2 p ε 0 E. ( 15
)
The physical predictions obtained by the GNOR and NHD models often differ substantially, as illustrated below. However, from a computational point of view the GNOR model only differs by the replacement β 2 → β 2 + D(γ -iω) in the frequency domain, whereby the nonlocal hydrodynamic parameter acquires an often non-negligible imaginary part. In the GNOR model we have the same additional boundary condition n • J = 0 on ∂Ω S as in the NHD model.

Specification to 2D TM mode

Now we can couple Maxwell's equation ( 1) with ( 9) for the NHD model, or similarly with [START_REF] Dressel | Verifying the Drude response[END_REF] for the GNOR model. We will compute light extinction by infinitely long nanowires. We take the wire axes along the z-direction and consider TM-polarized incident light, i.e. polarized in the (x, y)-plane. In this 2D setting, we can define E = (E x , E y ) T to be a vector and H = H z a scalar function. Coupling the time-harmonic Maxwell's equations and hydrodynamic Drude model ( 9), we have in 2D

         ∇ × H = -iωε 0 E + J, in Ω, ∇ × E = iωµ 0 H, in Ω, ∇(∇ • J) + ω(ω + iγ) β 2 J = iωω 2 p ε 0 β 2 E, in Ω S . (16)
If the Silver-Müller boundary condition (first-order absorbing boundary condition) [START_REF] Stupfel | Absorbing boundary conditions on arbitrary boundaries for the scalar and vector wave equations[END_REF] is applied on the boundary ∂Ω of the computational domain, then we have the boundary

conditions n × E -H = n × E inc -H inc = g inc , on ∂Ω, n • J = 0, on ∂Ω S , (17)
where E inc and H inc stand for the electromagnetic fields of the incoming light.

HDG formulations of nonlocal optical response models

The promise of hybridizable DG methods

In the Introduction some properties and advantages of DG and HDG methods were briefly mentioned, which we here explain in more detail. The classic DG method is seldomly employed for solving stationary problems, because it duplicates degrees of freedom (DOFs) on every internal edge. Thus the number of globally coupled DOFs is much greater than the number of DOFs required by conforming finite element methods for the same accuracy. Consequently, DG methods are expensive in terms of both CPU time and memory consumption. Hybridization of DG methods [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF] is devoted to addressing this issue while at the same time keeping all the advantages of DG methods. HDG methods introduce additional hybrid variables on the edges of the elements. Then we define the numerical traces arising from partial integration in the DG formulations through the hybrid variables. We can thus define the local (element-wise) solutions by hybrid variables. Conservativity conditions are imposed on numerical traces to ensure the continuity of the tangential component of the electric field and the normal component of the current density and to make the problem solvable. As a result, HDG methods produce a linear system in terms of the DOFs of the additional hybrid variables only. In this way, the number of globally coupled DOFs is greatly reduced as compared to the classic DG method. In a recent study [START_REF] Yakovlev | To CG or to HDG: A Comparative Study in 3D[END_REF], the authors showed that HDG methods outperform FEM in many cases.

Computational concepts and notations

In order to give a clear presentation of the HDG method, here we introduce some computational concepts and notations. We divide the computational domain Ω into triangle elements. The union of all the triangles is denoted by T h . By F h we denote the union of all edges of T h . Furthermore, F I h stands for the union of all the edges associated with the nanostructure. For an edge associated with two elements F = K + ∩ K -∈ F h , let (v ± , v ± ) be the traces of (v, v) on F from the interior of K ± , see Fig. 2, where we use the term trace to denote the restriction of a function on the boundaries of the elements [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF]. Note that from now on v is used to describe a general vector function instead of velocity. On every face, we define mean (average) values {•} and jumps • as

                   {v} F = 1 2 (v + + v -), {v} F = 1 2 (v + + v -), n × v F = n + × v + + n -× v -, n • v F = n + • v + + n -• v -, vt F = v + t + + v -t -,
where n ± denotes the outward unit norm vector to K ± and t ± denotes the unit tangent vectors to the boundaries ∂K ± such that t + × n + = 1 and t -× n -= 1. For the boundary

K + K - F v + v - n + n - Figure 2:
Two neighboring discretization elements (here: triangles) within the computational domain. An edge F is shared by two elements K + and K -. The outward normal vectors n + and n -point in opposite directions. A characteristic property of the DG method is that computed functions are allowed to be discontinuous across F (hence the "D" in DG). For example, for a function v, be it a scalar or a vector, its value on F from K + is v + , while its value on F from K -is v -, and these v + and v -are not necessarily equal. By contrast, the hybrid variables in the HDG method are single-valued on F .

edges, either on ∂Ω or on ∂Ω S , these expressions are modified as

               {v} F = v + , {v} F = v + , n × v F = n + × v + , n • v F = n + • v + , vt F = v + t + .
Let P p (D) denote the space of polynomial functions of degree at most p on a domain D. For any element K ∈ T h , let V p (K) be the space P p (K) and V p (K) the space (P p (K)) 2 . The discontinuous finite element spaces are then defined by

V p h = v ∈ L 2 (Ω) | v| K ∈ V p (K), ∀K ∈ T h , V p h = v ∈ (L 2 (Ω)) 2 | v| K ∈ V p (K), ∀K ∈ T h , (18)
where L 2 (Ω) is the space of square integrable functions on the domain Ω. We also introduce a traced finite element space

M p h = η ∈ L 2 (F h ) | η| F ∈ P p (F ), ∀F ∈ F h .
Note that M p h consists of functions which are continuous on an edge, but discontinuous at its ends. The restrictions of V p h , V p h and M p h in Ω S are denoted by V p h , V p h and M p h . For two vectorial functions u and v in (L 2 (D)) 2 , we introduce the inner product (u, v) D = D u•v dx, where • denotes the complex conjugation. Likewise for scalar functions u and v in L 2 (D), the inner product is defined as (u, v) D = D uv dx provided D is a domain in R 2 . Finally we define the edge overlap u, v F = F uv ds, where F is a specific edge. Accordingly, we can define the total edge overlap for the whole triangulation or for relevant subsets of edges. Important cases are

•, • F h = F ∈F h •, • F , •, • ∂Ω = F ∈F h ∩∂Ω •, • F , •, • F I h = F ∈F I h •, • F .
denoting, respectively, the total edge overlap on the computational domain, the cumulative edge overlap on the absorbing boundary of the computational domain, and finally the cumulative edge overlap on the nanostructure.

DG formulation of the coupled electrodynamical equations

We begin the construction of a DG implementation of the hydrodynamic Drude model by rewriting the coupled electrodynamical equations ( 16) into a system of first-order equations

             iωε 0 E + ∇ × H -J = 0 in Ω, iωµ 0 H -∇ × E = 0 in Ω, ∇q + γ -iω β 2 J - ω 2 p ε 0 β 2 E = 0 in Ω S , iωq -∇ • J = 0 in Ω S , (19)
where we introduced the scalar function q = (iω) -1 ∇ • J which coincides with a scaled charge density. In general, a DG method seeks an approximate solution (E h , H h , J h , q h ) in the space

V p h × V p h × V p h × V p
h that for each element K (in our case: for each discretization triangle) satisfies [START_REF] Bouajaji | High order discontinuous Galerkin method for the solution of 2D time-harmonic Maxwell's equations[END_REF]

               (iωε 0 E h , v) K + (∇ × H h , v) K -(J h , v) K = 0 ∀v ∈ V p (K), (iωµ 0 H h , v) K -(∇ × E h , v) K = 0 ∀v ∈ V p (K), (∇q h , v) K + γ -iω β 2 J h , v K - ω 2 p ε 0 β 2 E h , v K = 0 ∀v ∈ V p (K), (iωq h , v) K -(∇ • J h , v) K = 0 ∀v ∈ V p (K). (20)
The application of appropriate Green's formulas to this system of equations leads to terms on the element boundaries [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF]. These boundary terms are the keys to connect the elements, since the elements themselves are independent due to the nature of the discontinuous finite elements spaces of Eq. ( 18). In a DG method, one replaces the boundary terms by so-called numerical traces Êh , Ĥh , Ĵh and qh [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF][START_REF] Li | Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell's equations[END_REF], which are also known as 'numerical fluxes' in the literature [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF]. These numerical traces are defined as

           Ĥh = {H h } + α E n × E h , n × Êh = {n × E h } + α H H h , qh = {q h } + α J n • J h , n • Ĵh = {n • J} + α q q h . (21)
In these definitions there is still freedom to choose values for the α parameters, and this corresponds to different DG schemes: by setting α E = α H = α J = α q = 0, one obtains the centered flux DG scheme. With α E = α H = α J = α q = 1, one obtains the upwind flux DG scheme [START_REF] Bouajaji | High order discontinuous Galerkin method for the solution of 2D time-harmonic Maxwell's equations[END_REF]. For more validated DG schemes, we refer the interested readers to Ref. [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF].

Having defined the numerical traces, we finally form a global system of linear equations involving all the DOFs on all the elements

               (iωε 0 E h , v) K + (H h , ∇ × v) K -Ĥh , n × v ∂K -(J h , v) K = 0 ∀v ∈ V p (K), (iωµ 0 H h , v) K -(E h , ∇ × v) K -n × Êh , v ∂K = 0 ∀v ∈ V p (K), -(q h , ∇ • v) K + qh , n • v ∂K + γ -iω β 2 J h , v K - ω 2 p ε 0 β 2 E h , v K = 0 ∀v ∈ V p (K), (iωq h , v) K + (J h , ∇v) K -n • Ĵh , v ∂K = 0 ∀v ∈ V p (K), (22)
which are coupled equations that are valid whatever DG scheme is adopted.

Hybridizable DG implementation of the electrodynamical equations

In Sec. 3.1 we mentioned that hybridized DG methods have advantages as compared to the classic DG schemes, and here we discuss the hybridized approach in more detail. Unlike in the above classic DG formulations where the numerical traces directly couple the values from the elements on both sides of the edges, in a HDG formulation the numerical traces are defined through hybrid variables. Introducing two hybrid variables λ h and η h which live only on the boundaries of the elements, we define the numerical traces by

           Ĥh = λ h , Êh = E h + τ λ (λ h -H h )t, qh = η h , Ĵh = J + τ η (q h -η h )n, (23)
where τ λ and τ η are two stabilization parameters. Replacing the numerical traces in [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF] with the expressions in [START_REF] Li | A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell's equation[END_REF] and applying Green's formulas to the first and fourth equations in [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF], we obtain the local formulation of the HDG method as

               (iωε 0 E h , v) K + (H h , ∇ × v) K -λ h , n × v ∂K -(J h , v) K = 0, ∀v ∈ V p (K), (iωµ 0 H h , v) K -(∇ × E h , v) K + τ λ (H h -λ h ), v ∂K = 0, ∀v ∈ V p (K), -(q h , ∇ • v) K + η h , n • v ∂K + γ -iω β 2 J h , v K - ω 2 p ε 0 β 2 E h , v K = 0, ∀v ∈ V p (K), (iωq h , v) K -(∇ • J h , v) K -τ η (q h -η h ), v ∂K = 0, ∀v ∈ V p (K). ( 24
)
One can solve the local fields element by element once the solutions for λ h and η h are obtained. In order to make the problem solvable, we need to employ global conditions

n × Êh , v F h -λ h , v ∂Ω = g inc,v ∂Ω , ∀v ∈ M p h , n • Ĵh , v F I h = 0, ∀v ∈ M p h . (25)
The first relation in [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF] weakly enforces the continuity condition for the tangential component of the electric field across any edges, and also takes into account the Silver-Müller absorbing boundary condition. The other global condition in Eq. ( 25) weakly enforces the continuity condition for the normal component of the current density across any edges. The additional boundary condition on the surface of the nanostructure is implicitly contained in this relation. Substituting Êh and Ĵh in [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF] with the definitions in [START_REF] Li | A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell's equation[END_REF], we arrive at the global reduced system of equations

n × E h -τ λ (λ h -H h ), v F h -λ h , v ∂Ω = g inc , v ∂Ω , ∀v ∈ M p h , n • J h + τ η (q h -η h ), v F I h = 0, ∀v ∈ M p h . ( 26
)
Note that we used the fact that n × t = -1 in [START_REF] Huang | Theoretical and numerical analysis of a non-local dispersion model for light interaction with metallic nanostructures[END_REF]. The two relations in Eq. ( 26) are not independent. They are coupled through the local solutions of E, H, J and q of the local equations ( 24). Remark I. The proposed HDG formulation for the global system ( 26) is naturally consistent with the boundary conditions, both on the artificial boundary and on the medium boundary.

Remark II. Globally, we only need to solve Eq. ( 26), in which the fields E h , H h , J h and q h are replaced by the solutions in terms of λ h and η h from the local problems [START_REF] He | Optimized Schwarz algorithms for solving time-harmonic Maxwell's equations discretized by a hybridizable discontinuous Galerkin method[END_REF]. So the global DOFs are associated with λ h in the whole computational domain, while they are associate with η h only within in the material medium. The discretization leads to a system of linear equations

A λ h η h = g inc h 0 , ( 27
)
where λ h and η h are vectors accounting for the degrees of freedom of the hybrid variables λ and η respectively, and the coefficient matrix A is large and sparse.

Numerical results

In this section we present numerical results to validate the proposed HDG formulations. All HDG methods have been implemented in Fortran 90. All our tests are performed on a Macbook with a 1.3 GHz Inter Core i5 CPU and 4 GB memory. We employ the multifrontal sparse direct solver MUMPS [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] to solve the discretized systems of linear equations.

In HDG methods, we calculate the total fields E tot and H tot . The scattered fields are then calculated by subtracting the incident field from the total fields. We use HDG-P p to denote the HDG method with interpolation order p. Here we choose fixed values τ λ = τ η = 1 for the stabilization parameters. Different choices are discussed in Ref. [START_REF] Gopalakrishnan | Stabilization in relation to wavenumber in HDG methods[END_REF].

Convergence study: Wave propagation in a cavity

While elsewhere in this article we focus on nanowire structures, here we first study the convergence of our method by considering wave propagation in a cavity. This cavity is assumed to be a square domain Ω = {(x, y) ∈ [0, L] × [0, L]} with the PEC boundary condition and hard-wall condition n × E = 0, and n • J = 0, on ∂Ω .

This test case can be viewed as the frequency-domain version of the first test case in [START_REF] Schmitt | A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects[END_REF]. The simplicity is achieved by introducing artificial current density and electric field, such that the analytical solutions coincide with Maxwell's equations and with the hydrodynamic equation

             ikE + ∇ × H = J -J a , ikH -∇ × E = 0, ∇q - iω β 2 J = - γ β 2 J + γ β 2 J a + ω 2 p ε 0 β 2 E - ω 2 p ε 0 β 2 E a , iωq -∇ • J = 0, (28)
where k = ω c is the wave number, with c being the light speed. We make this modification to unify the scale of the electric and magnetic fields. The artificial terms J a and E a are also the analytical solution to this equation ( 28):

E a = √ 2 2 i -cos( √ 2 2 kx) sin( √ 2 2 ky) sin( √ 2 2 kx) cos( √ 2 2 ky)
,

J a = - √ 2µ 0 kβ 2 2ω sin( √ 2 2 kx) cos( √ 2 2 ky) cos( √ 2 2 kx) sin( √ 2 2 ky) . (29)
We only take the real part of H and J and the imaginary part of E and q into consideration.

In order to have this analytical solution [START_REF] Raza | Nonlocal study of ultimate plasmon hybridization[END_REF], one needs to set the length of the square

L = √ 2π k and β 2 = ω 2 k 2 .
The convergence history of the HDG method with interpolation order P p (p = 1, 2, 3) is given in Table 1 and Figure 3. Mesh size h is the edge length of elements associated to the boundary ∂Ω . The convergence orders are calculated by

log( E a -E h h 2 Ω / E a -E h h 1 Ω ) log(h 2 /h 1 )
,

where h 1 and h 2 denote a coarse and a refined mesh size, respectively. From Table 1 and Figure 3, we observe that the proposed HDG method has an optimal convergence order which is p + 1 for HDG-P p .

Benchmark problem: a cylindrical plasmonic nanowire

As our benchmark problem we consider the plasmonic behavior of a cylindrical nanowire. This has been used as a convenient benchmark problem for other numerical methods before [START_REF] Toscano | Modified field enhancement and extinction by plasmonic nanowire dimers due to nonlocal response[END_REF][START_REF] Hiremath | Numerical solution of nonlocal hydrodynamic Drude model for arbitrary shaped nano-plasmonic structures using Nédélec finite elements[END_REF] because analytical solutions exist both for the local and for the NHD models, see the derivation in Ref. [START_REF] Ruppin | Extinction properties of thin metallic nanowires[END_REF]. We make use of the fact that the analytical Mie solution of Ref. [START_REF] Ruppin | Extinction properties of thin metallic nanowires[END_REF] allows making the nonlocal parameter β complex-valued. This enables us to benchmark our HDG simulations against exact analytical results for the GNOR model as well. (For comparison, optical properties of a sphere in the GNOR model, based on exact Mie results, are discussed in Ref. [START_REF] Raza | Nonlocal optical response in metallic nanostructures[END_REF].) For the NHD model, the configuration of the nanowire is taken to be the same as that in the first test in [START_REF] Hiremath | Numerical solution of nonlocal hydrodynamic Drude model for arbitrary shaped nano-plasmonic structures using Nédélec finite elements[END_REF]: the radius of the cylinder is 2 nm, no interband transitions are considered, the plasma frequency ω p = 8.65 × 10 15 , the damping constant γ = 0.01ω p , the Fermi velocity v F = 1.07 × 10 6 , and β 2 = 3 5 v 2 F . For the GNOR model, we use the same parameters and furthermore we take D = 2.04 × 10 -4 [START_REF] Mortensen | A generalized non-local optical response theory for plasmonic nanostructures[END_REF]. An artificial absorbing boundary is set to be a concentric circle with a radius of 100 nm.

As our benchmark observable we will calculate the Extinction Cross Section (ECS, σ ext ), which is given by the sum of the scattering cross section σ sca and the absorption cross section σ abs [START_REF] Berg | General derivation of the total electromagnetic cross sections for an arbitrary particle[END_REF],

σ ext = σ sca + σ abs .
More precisely, for the cylindrical nanowire we consider the extinction cross section per wire length, which actually has the units of a length. We scale this quantity by the diameter 2r of the nanowire to obtain a dimensionless normalized extinction cross section that we denote by σ ext . It can be expressed as the sum of scaled scattering and absorption cross sections,

σ sca = 1 2r Re S (E sca × H sca ) • n dS, and σ abs = - 1 2r Re S (E tot × H tot ) • n dS.
Here the integrations are performed along a closed path around the nanowire, and Re denotes the real part.

The simulations are performed on a mesh with 4,513 nodes, 8,896 elements and 13,280 edges of which 722 edges are located inside the nanostructure. The ECS is presented in Figure 4. Curvilinear treatment is employed for high-order accuracy, where the curved edges are geometrically approximated by second-order curves instead of straight lines [START_REF] Li | Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell's equations[END_REF]. From Figure 4 we can observe that the fourth-order HDG method produces an ECS curve that matches the analytical solution very well. By contrast, the first-order method is not accurate enough on this mesh. Contour plots of the electric field and the current density are presented in Figure 5. These results match well with corresponding results in Ref. [START_REF] Hiremath | Numerical solution of nonlocal hydrodynamic Drude model for arbitrary shaped nano-plasmonic structures using Nédélec finite elements[END_REF] despite the lower resolution, probably because our simulation is performed on a coarser mesh. Comparing the two subfigures in Figure 4, we also find that the ECS curve for the GNOR model is smoother than for the NHD model. But this has a physical rather than a numerical origin. In particular the standing bulk plasmon resonances above the plasma frequency in the NHD model are essentially washed out by the introduced diffusion in the GNOR model. The ECS curves of HDG-P 2 and HDG-P 3 are are not presented in Figure 4, but we found that they they lie in between the displayed curves of HDG-P 1 and HDG-P 4 .

In our 2D simulations, we use a sparse direct solver MUMPS [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] to solve the resulting systems of linear equations. We need to solve a linear system at each frequency. The computational performance mainly relies on the size of the coefficient matrices, i.e. the number of degrees of freedom (#DOF). The computational performance for one frequency 2 we can see that the HDG-P 4 is more expensive than HDG-P 1 in CPU time for both construction and factorization. However, high-order methods are preferable because they costs less for the same accuracy [START_REF] Li | Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell's equations[END_REF].

Dimer of cylindrical nanowires

Plasmonic dimer structures with small gaps are both experimentally interesting and computationally challenging because of high field enhancements in the gap region [START_REF] Toscano | Modified field enhancement and extinction by plasmonic nanowire dimers due to nonlocal response[END_REF][START_REF] Gallinet | Numerical methods for nanophotonics: standard problems and future challenges Laser[END_REF][START_REF] Raza | Nonlocal optical response in metallic nanostructures[END_REF]. Here we present our HDG simulations of a cylindrical gold dimer geometry as shown in Figure 6(a), and this particular configuration is from Ref. [START_REF] Raza | Nonlocal study of ultimate plasmon hybridization[END_REF]. A typical mesh is shown in Figure 6(b). On a mesh with 5,829 nodes, 11,520 triangles and 17,348 edges with 3,712 edges inside the nanostructure, we calculate the ECS curve by HDG-P 4 . The size of matrix for HDG-P 1 is 105, 300 × 105, 300, the matrix construction CPU time is 5.2 seconds, the factorization CPU time is 6.9 seconds for one frequency, and the memory cost is 717 MB. For the material properties gold we use the same values as in Ref. [START_REF] Toscano | Modified field enhancement and extinction by plasmonic nanowire dimers due to nonlocal response[END_REF]: the plasma frequency ω p = 1.34 × 10 16 , damping constant γ = 1.14 × 10 14 , the Fermi velocity v F = 1.39 × 10 6 , and the nonlocal parameter β is determined by β 2 = 3 5 v 2 F . The incoming plane wave of light is incident perpendicular to the line connecting the centers of the two circles, with a linear polarization parallel to this line (TM-or p-polarization). A comparison of the ECS curves is presented in Figure 7. Overall, there are small but clear differences, illustrating that nonlocal response effects occur even for dimer structures for which the corresponding monomers (r = 30 nm nanowires) would show essentially no nonlocal effects [START_REF] Toscano | Modified field enhancement and extinction by plasmonic nanowire dimers due to nonlocal response[END_REF]. Both nonlocal models have blueshifted resonances as compared to the local model, and resonances in the GNOR model are less pronounced than in the local and NHD models. For smaller gap sizes, nonlocal blueshifts are larger and resonances are broadened more (the latter only in the GNOR model). Field distributions at the same particular frequency for the NHD and the GNOR models are compared in Figure 8. The figure illustrates the generic features that the GNOR model washes out some finer details of the field distributions, and also that minimal and maximal field values lie closer together in the GNOR model.

Conclusions

This paper introduces a HDG method to solve the nonlocal hydrodynamic Drude model and the GNOR model, both of which are often employed to describe light-matter interactions of nanostructures. The numerical fluxes are expressed in terms of two newly introduced hybrid terms. Only the hybrid unknowns are involved in the global problem. The local problems are solved element-by-element once the hybrid terms are obtained. The proposed HDG formulations naturally couple the hard-wall boundary condition. Numerical results indicate that the HDG method converges at the optimal rate. Our benchmark simulations for a cylindrical nanowire and our calculations for a dimer structure show that the HDG method is a promising method in nanophotonics. Building on these results, in the near future we plan to generalize our computations to 3D structures, and to introduce domain decomposition and model order reduction into nanophotonic computations.

Figure 3 :

 3 Figure 3: Convergence history of the proposed HDG method for the cavity problem.

Figure 4 :

 4 Figure 4: Extinction cross section of a Na cylinder with 2 nanometer radius in a free-space background, for a TM-polarized normally incident plane wave. The cylinder is described both by (a) the NHD model and by (b) the GNOR model. The simple wire geometry serves as an excellent benchmark problem: Analytically exact calculations (nonlocal Mie theory) are compared with HDG methods of different interpolation order.

Figure 5 :

 5 Figure 5: The electric-field and current-density distributions of the light-matter interaction of a Na nanowire.In the upper row, we show the distributions on the fourth-order nonlocal resonance at ω/ω p = 1.227 for the NHD model. For comparison, in the bottom row we show the corresponding distributions for the GNOR model.

Figure 6 :

 6 Figure 6: A cylindrical gold dimer nanowire. Panel (a): The geometry of the nanowire dimer. An imaginary rectangle is introduced around the dimer for the calculation of the cross section. Both nanowires have a radius r = 30 nm and their gap distance is d = 3 nm. Panel (b): A typical mesh. The large circle is the artificial absorbing boundary.

Figure 7 :

 7 Figure 7: Comparison of extinction cross sections of a gold dimer as calculated with the local Drude model, the NHD model and the GNOR model. The configuration is shown in Figure 6 and the material parameters are given in subsection 4.3.

Figure 8 :

 8 Figure 8: Various field distributions in the gold dimer when illuminated by a plane wave of light. On the top line we show the distributions at the third SPR of the NHD model, at ω/ω p = 0.66. On the second line we show the corresponding distributions in the GNOR model at the same frequency.

Table 1 :

 1 Convergence results for the cavity problem.

		HDG-P 1		HDG-P 2		HDG-P 3
	h	error	order	error	order	error	order
	5 -8	1.67 × 10 -9	-	4.52 × 10 -10	-	2.04 × 10 -11	-
	2.5 -8 4.10 × 10 -10

Table 2 :

 2 Computational performance of the nanowire problem.

	#DOF t construction (second) t factorization (second) memory (MB)
	HDG-P 1 28,260	0.067	0.36	74
	HDG-P 4 70,650	2.4	3.3	418
	is given in

Table 2 ,

 2 where t construction denotes the CPU time for construction the matrices, t factorization denotes the CPU time used by MUMPS for the factorization of the coefficient matrix A (27), and memory denotes the memory consumed by MUMPS. From Table

https://www.lumerical.com/

http://www.comsol.com/

Acknowledgments

The first author was supported by the NSFC (11301057) and the Fundamental Research Funds for the Central Universities (ZYGX2014J082). N. A. M. and M. W. acknowledge support from the Danish Council for Independent Research (FNU 1323-00087). M. W. acknowledges support from the Villum Foundation via the VKR Centre of Excellence NATEC-II. The Center for Nanostructured Graphene is sponsored by the Danish National Research Foundation, Project DNRF103.


H Piret

P Granjon

N Guillet

V Cattin

Tracking of electrochemical impedance of batteries

Keywords: battery monitoring system, electrochemical impedance spectroscopy, frequency domain estimation, wideband identification, recursive estimation algorithm

This paper presents an evolutionary battery impedance estimation method, which can be easily embedded in vehicles or nomad devices. The proposed method not only allows an accurate frequency impedance estimation, but also a tracking of its temporal evolution contrary to classical electrochemical impedance spectroscopy methods. Taking into account constraints of cost and complexity, we propose to use the existing electronics of current control to perform a frequency evolutionary estimation of the electrochemical impedance. The developed method uses a simple wideband input signal, and relies on a recursive local average of Fourier transforms. The averaging is controlled by a single parameter, managing a trade-off between tracking and estimation performance. This normalized parameter allows to correctly adapt the behavior of the proposed estimator to the variations of the impedance. The advantage of the proposed method is twofold: the method is easy to embed into a simple electronic circuit, and the battery impedance estimator is evolutionary. The ability of the method to monitor the impedance over time is demonstrated on a simulator, and on a real Lithium ion battery, on which a repeatability study is carried out. The experiments reveal good tracking results, and estimation performance as accurate as the usual laboratory approaches.

Introduction

The recent and future expansion of electric vehicles or nomad devices inevitably leads to the development of efficient battery management systems (BMS). Such systems must continuously determine the state of the monitored battery from several measurements. However, in order to preserve the battery integrity, only non invasive and non destructive measurement methods are used, and most BMS measure external quantities such as the current flowing through the battery, the voltage across its terminals and its surface temperature [START_REF] Coleman | State-of-charge determination from emf voltage estimation: Using impedance, terminal voltage, and current for lead-acid and lithium-ion batteries[END_REF] [2]. From these measurements, one way to obtain interesting information representative of the current state of the battery is to estimate its electrochemical impedance [START_REF] Howey | Impedance measurement for advanced battery management system[END_REF]. Indeed this quantity describes the dynamic behavior of the battery and regularly changes with the evolution of its internal temperature, state of charge (SoC) and state of health (SoH) [START_REF] Waag | Experimental investigation of the lithiumion battery impedance characteristic at various conditions and aging states and its influence on the application[END_REF]. Consequently, the electrochemical impedance is used in many methods to estimate the internal temperature [START_REF] Zhu | A new lithuim-ion battery internal temperature in-line estimate method based on electrochemical impedance spectroscopy measurement[END_REF] [START_REF] Richardson | Battery internal temperature estimation by combined impedance and surface temperature measurement[END_REF] [7], the SoC [START_REF] Bundy | An electrochemical impedance spectroscopy method for prediction of the state of charge of a nickel metal battery at open circuit and during discharges[END_REF] [9] [START_REF] Blanke | Impedance measurements on leadacid batteries for state-of-charge, state-of-health and cranking capability prognosis in electric and hybrid electric vehicles[END_REF] and the SoH [START_REF] Blanke | Impedance measurements on leadacid batteries for state-of-charge, state-of-health and cranking capability prognosis in electric and hybrid electric vehicles[END_REF] [START_REF] Troltzsch | Characterizing aging effects of lithium ion batteries by impedance spectroscopy[END_REF] of the monitored battery. Usual procedures used to estimate a battery electrochemical impedance belong to the class of active identification methods [START_REF] Ljung | System Identification -Theory for the User[END_REF]: known variations are added to the battery input, and the corresponding output variations are measured and used to estimate the unknown impedance on a given frequency band. In the rest of this paper, the battery impedance is estimated in galvanostatic mode: the current i(t) flowing through the battery is considered as the input and the voltage u(t) across the battery terminals as the output. This choice can be justified by the fact that the current can be easily driven by very simple and cheap electronic devices (such as a transistor for example), which is a strong requirement for embedded systems for which this work is developed. Active identification methods rely on two main assumptions:

• Firstly, variations of the additional current used to estimate the impedance are chosen sufficiently small for the battery to have a linear behavior with respect to these variations. Under this assumption, the battery can be considered as a linear system.

• Secondly, parameters on which the battery characteristics depend are assumed to remain constant during the measurement process. Under this assumption, the battery can be considered as a time-invariant system during the measurement time.

Jointly, these two assumptions allow to consider the battery as a linear and time-invariant (LTI) system regarding the additive input-output variations and during the measurement time. In that case, the battery admits a well defined frequency response function, corresponding to its electrochemical impedance Z(f ) and verifying the following frequency relationship:

Z(f ) = S ui (f ) S ii (f ) if S ii (f ) = 0. (1)
In this equation, S ui (f ) is the cross power spectral density (CPSD) between voltage and current variations, while S ii (f ) is the power spectral density (PSD) of current variations only [START_REF] Shin | Fundamentals of Signal Processing for Sound and Vibration Engineers[END_REF][START_REF] Pintelon | System Identification -a Frequency Domain Approach[END_REF]. Eq. ( 1) highlights the fact that Z(f ) can only be estimated in the frequency bands where the input current variations contain power, i.e. where their PSD is different from 0.

The validity of the LTI assumption for the battery and equivalently the validity of Eq. ( 1) can be checked by using the notion of magnitude squared spectral coherence [START_REF] Bendat | Random Data -Analysis and Measurement Procedures[END_REF] defined as:

C ui (f ) = |S ui (f )| 2 S uu (f )S ii (f ) , (2)
where S uu (f ) is the PSD of voltage variations.

This frequency domain function is a statistical quantity normalized between 0 and 1, that can be interpreted as the magnitude squared correlation coefficient between the spectral components of the voltage and the current around a given frequency f . It gives a normalized measurement of how linearly the spectral components of these two signals are related to each other. It has been shown

for example in [START_REF] Bendat | Random Data -Analysis and Measurement Procedures[END_REF] that in case of low measurement noise, the identified system can be considered as LTI in frequency bands where the magnitude squared coherence is close to 1, while the LTI assumption can be rejected in frequency bands where it remains close to 0. This quantity has been used to check the 3 LTI assumption for batteries in [START_REF] Al-Nazer | Broadband identification of battery electrical impedance for hevs[END_REF].

Several methods have been developed to estimate the impedance, such as the time domain step response. The disadvantage of this approach is related to the measured response to an impulse or step input that often have a small amplitude in comparison to the noise, especially when the battery impedance is low. That is why this technique requires extra-large inputs to reach good signal to noise ratios and finally good estimation performance. These large inputs often induce non-linear behavior, explaining why this approach is not used in this study.

One of the authoritative methods for battery impedance measurements is the narrowband electrochemical impedance spectroscopy (FFT-EIS) [17][18]. In this method, a single sine wave with low amplitude and fixed frequency is used as input signal. Eq. ( 1) is then valid at the sine frequency only, and Z(f ) can be estimated at that particular frequency only, justifying the term "narrowband".

If Z(f ) must be estimated for several frequencies, the same measurement process has to be sequentially done for each desired frequency. An efficient way to avoid this sequential implementation and estimate the impedance for a discrete set of frequencies at one time is to use a multisine approach [START_REF] Pintelon | System Identification -a Frequency Domain Approach[END_REF]. In that case, the input signal consists of a sum of sines which frequencies correspond to the desired set.

For wideband method, input signals are wideband in the sense that their PSD is different from 0 on a continuous frequency band. In that case, Eq. ( 1) is valid all over that frequency band, where Z(f ) can be estimated whatever f . Several options are available to choose a wideband input signal for system identification, the most popular being swept sines, random noises and pseudorandom binary sequences (PRBS) [START_REF] Pintelon | System Identification -a Frequency Domain Approach[END_REF].

Eq. ( 1) and ( 2) clearly show that the estimation of Z(f ) relies exclusively on basic spectral quantities such as PSD and CPSD. A simple and efficient estimator usually used for such quantities is the Welch modified periodogram [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF]. The signals are first divided into L consecutive blocks of same length by using a time window. The discrete Fourier transform (DFT) of each block of data is then computed by using a fast Fourier transform algorithm. Finally, the L obtained result. As an example, Eq. ( 3) gives the expression of the voltage-current CPSD estimator Ŝui (f ):

Pui k (f ) = AU k (f )I k * (f ), (3a) Ŝui (f ) = 1 L L-1 k=0 Pui k (f ), ( 3b
)
where A is a normalization factor, * denotes complex conjugation, and U k (f ) (I k (f ) respectively) is the DFT of the k th block of voltage (current respectively) signal.

In Eq. (3a), Pui k (f ) is the cross-periodogram of the k th blocks of voltage and current signals, and Eq. (3b) clearly shows that the estimated CPSD is given by an arithmetic averaging of the L cross-periodograms obtained from the acquired data. Obviously, same type of estimators can be obtained for the current and voltage PSDs Ŝii (f ) and Ŝuu (f ) by using exclusively U k (f ) or

I k (f ) in Eq. (3).
A simple impedance estimator Ẑ(f ) is obtained by using Eq. (3) in the impedance definition given by Eq. [START_REF] Coleman | State-of-charge determination from emf voltage estimation: Using impedance, terminal voltage, and current for lead-acid and lithium-ion batteries[END_REF].

Ẑ(f ) = Ŝui (f ) Ŝii (f ) where Ŝii (f ) = 0. ( 4
)
Following the same principle, Eq. (3) used in Eq. ( 2) leads to the estimator of the magnitude squared spectral coherence between battery voltage and current:

Ĉui (f ) = | Ŝui (f )| 2 Ŝuu (f ) Ŝii (f ) . (5)
Therefore, Eq. ( 3), ( 4) and ( 5) form together the battery identification algorithm:

• Eq. ( 3) and ( 4) give access to the battery impedance estimate in the frequency band of interest,

• Eq. ( 3) and ( 5) can be used to check the validity of the LTI assumption for the battery behavior, or equivalently the validity of the estimated impedance.

The performance of this identification strategy has been studied in numerous general books such as [START_REF] Ljung | System Identification -Theory for the User[END_REF], [START_REF] Pintelon | System Identification -a Frequency Domain Approach[END_REF] and [START_REF] Bendat | Random Data -Analysis and Measurement Procedures[END_REF], and in more specific works such as [START_REF] Al-Nazer | Broadband identification of battery electrical impedance for hevs[END_REF] [20] for battery impedance estimation. It has been shown that on the one hand, the variance of the impedance estimator is inversely proportional to the number L of averaged periodograms in Eq (3), which is a quite intuitive result. On the other hand, this variance is also proportional to 1-Cui(f ) Cui(f ) showing that the best estimation performance will be reached when C ui (f ) is close to 1, that is for true LTI systems and small measurement noises.

While single sine and multisine FFT-EIS reach very accurate estimation performance [START_REF] Pintelon | System Identification -a Frequency Domain Approach[END_REF], the set of frequencies at which Z(f ) is estimated depends on the chosen input signal. In other words, each time Z(f ) must be estimated at slightly different frequencies, a new input signal must be designed and the whole measurement process must be repeated. Moreover, the generation of sine waves with different frequency or multisine signals requires advanced electronic generators. These disadvantages concern both EIS based on phase sensitive detector [START_REF] Lasia | Electrochemical Impedance Spectroscopy and its Applications[END_REF] and FFT-EIS. Even if the generation of input signals is much easier for square current EIS [START_REF] Yokoshima | Application of electrochemical impedance spectroscopy to ferri/ferrocyanide redox couple and lithium ion battery systems using a square wave as signal input[END_REF], this method only allows to estimate the impedance for a small discrete set of frequency: the fundamental frequency of the square wave and its odd harmonics.

Wideband methods relying on PRBS signals are more suitable for embedded systems. Indeed, such two-level signals can be generated with very simple electronic circuits. PRBS signals reach good estimation performance for electrochemical impedance of batteries [START_REF] Al-Nazer | Broadband identification of battery electrical impedance for hevs[END_REF], although it leads to slightly less accurate results than multisine signals [START_REF] Pintelon | System Identification -a Frequency Domain Approach[END_REF]. Moreover, the same measurement data set can be used to estimate Z(f ) whatever f in the frequency band of the chosen input signal.

However, a joint property of these identification strategies limits their application to embedded systems: they only lead to one estimate of Z(f ) per measurement. Consequently, each time a new estimate of Z(f ) is needed, the whole measurement process must be repeated, strongly limiting the capacity to follow the time evolution of Z(f ) that can be significant in case of strongly varying operating conditions.

The method detailed in the next section is able to efficiently track the variations of Z(f ) over time thanks to a recursive implementation of the wideband approach.

In this context, this work aims at developing a method dedicated to battery electrochemical impedance estimation and tracking, which is sufficiently efficient to precisely track the temporal variations of this quantity, and suitable to be embedded in a vehicle or a nomad device. The proposed method and the corresponding results are presented in this paper with the following organization.

Section 2 presents the proposed algorithm based on a local averaging strategy recursively implemented in the frequency domain, and the influence of its main parameters. This algorithm is tested on a simulator and on an experiment conducted on a real battery. The section 3 describes the protocol used in each case.

The section 4 details and discusses both the simulate and experimental results.

General principle

Time-varying impedance spectroscopy

The goal of this study is to obtain a time-dependent estimator of the battery electrochemical impedance. This estimator must be sufficiently accurate to give exact estimates of Z(f ), but must also have tracking capabilities to follow its eventual variations over time. One possibility is to use the wideband principles previously detailed, and adapt them to the time-varying case. In the wideband approach, the impedance estimation entirely relies on the Welch modified periodogram summarized in Eq. (3). More particularly, Eq. (3b) shows that this spectral estimator is based on a global average of the whole set of pediodograms obtained from the acquired data. The main idea behind the proposed method is to replace this global average with a local average, leading to time-dependent estimators for the different spectral quantities, and consequently for the impedance. In the next paragraph, two recursive local averaging strategies are detailed and compared.

Local averaging strategies

Among all averaging strategies, the ones that can be implemented recursively are particularly suitable for embedded systems. In this section, only the two simplest ones are detailed and compared: sliding window averaging, and exponential averaging.

The time-varying spectral estimator obtained with a sliding window averaging consists in the average of the N last periodograms. Following notations used in Eq. ( 3), the corresponding estimator for the voltage-current CPSD around the block of data number k is given by:

Ŝui k (f ) = 1 N k n=k-N +1 Puin (f ), (6a)
= Ŝui k-1 (f ) + Pui k (f ) -Pui k-N (f ) N , (6b)
where the number of averaged periodograms N is smallest than the total number of periodograms L. Eq. (6b) corresponds to the recursive implementation of this estimator, which is more suitable for embedded systems. Obviously, the averaging effect increases with the size of the sliding window N .

The spectral estimator obtained with an exponential averaging is given by the following relation:

Ŝui k (f ) = α Ŝui k-1 (f ) + (1 -α) Pui k (f ), (7)
where α ∈ [0, 1[ is usually called the forgetting factor. Eq. ( 7) clearly shows that when α = 0, the voltage-current CPSD estimate for block k Ŝui k (f ) is simply the cross-periodogram of the corresponding data blocks Pui k (f ). In that case, no average is performed, leading to poor estimation performance regarding measurement noise. On the contrary, the estimated CPSD is more and more averaged when α increases and tends to 1.

As noticed for Eq. ( 3), similar time-varying estimates are obtained for PSDs by replacing cross-periodograms with periodograms in Eq. ( 6) and [START_REF] Schmid | Measurement of the internal cell temperature via impedance: Evaluation and application of a new method[END_REF].

Two different criteria are used to compare the previous strategies: their estimation performance, and their implementation complexity.

Estimation performance can be easily studied by considering each strategy as a digital filter. Indeed, Eq. (6b) corresponds to the recursive equation of a causal finite impulse response (FIR) filter of order N with all weights equal to 1 N , while Eq. ( 7) is the recursive equation of a causal infinite impulse response (IIR) filter of order 1 with one pole α. Then, the following filter properties can be used to characterize the performance of each strategy:

• Response time at r%: this quantity is the number of samples from which the step response of the considered filter reaches r% of it's final value. In the present context, it can be interpreted as the convergence time at r% of the averaging strategy.

• Equivalent noise bandwidth: this quantity corresponds to the variance (i.e. the amount of residual fluctuations) at the filter output when supplied with a white noise of variance 1. For lowpass filters as in this study, it represents the significance of the averaging effect realized by the filter: strong averaging effect for values close to 0, and weak averaging effect for values close to one.

The values of these characteristics are determined for each averaging strategy by using the weights of the corresponding filter, and are given in Table . 1. Obviously, these characteristics exclusively depend on parameters N and α, and show that the two strategies can reach similar estimation performance depending on the values of these two independent parameters. These results are consequently not useful to determine the best strategy, but give interesting links between N and α. For example, α can be set such that Eq. ( 7) has the same averaging effect as Eq. ( 6). In that case, the two averaging strategies must have the same equivalent noise bandwidth, leading to the following equality:

α = N -1 N + 1 . (8)
The algorithm complexity can be directly deduced from the recursive equations (6b) and [START_REF] Schmid | Measurement of the internal cell temperature via impedance: Evaluation and application of a new method[END_REF]. The implementation of these two equations leads to a similar number of operations, but the sliding window approach clearly needs a larger amount of memory. Indeed, in Eq. ( 7) only two values must be stored to update the spectral estimator whatever the value of α. On the contrary, the memory depth needed for Eq. (6b) is N + 2, leading to a large amount of memory for large N . Therefore, the exponential averaging is much more interesting for embedded systems regarding this second criterion, which justifies the choice of this strategy in the final estimator.

Time-varying impedance estimator

As previously mentioned, the global structure of the proposed algorithm is the same as for the wideband approach, and is described in Fig. 1. The main difference concerns the estimation of spectral quantities which are now locally averaged around each block of data with an exponential averaging strategy. The final expression of the impedance estimated for the block number k is then:

Ẑk (f ) = Ŝui k (f ) Ŝii k (f ) , (9)
where the spectral quantities Ŝui k (f ) and Ŝii k (f ) are updated thanks to Eq. ( 7).

The obtained estimator is time-varying and updated at each new block of data in order to correctly estimate the true value of the battery impedance, and efficiently track its eventual variations over time. On the one hand, its tracking performance is related to the convergence time or similarly the response time of the averaging strategy: the smaller the convergence time is, the higher the tracking performance is. On the other hand, the estimation performance is directly related to the amount of residual variations after averaging, and therefore to the equivalent noise bandwidth of the averaging strategy: the lower the equivalent noise bandwidth is, the higher the estimation performance is. From expressions given in Table . 1, it is clear that the tradeoff between these two types of performance is controlled by the parameter α:

• When α is small and close to 0, the convergence time is small leading to high tracking capabilities. But in this case, the equivalent noise bandwidth is high meaning poor averaging effect and poor estimation performance.

• On the contrary if α is high and close to 1, the convergence time is large and the tracking capabilities decrease. However, the equivalent noise bandwidth is small in that case, leading to high averaging effect and good estimation performance.

Finally, Eq. ( 8) gives an interesting way to set the value of α: if the user wants the exponential averaging be equivalent to a sliding window of size N , α must be set by using this relation.

Following the same idea as for impedance estimation, the global estimator of the magnitude squared spectral coherence given in Eq. ( 5) can be adapted to the time-varying case by using local exponential averaging estimates of PSD and CPSD as shown in Eq. [START_REF] Blanke | Impedance measurements on leadacid batteries for state-of-charge, state-of-health and cranking capability prognosis in electric and hybrid electric vehicles[END_REF].

Ĉui k (f ) = | Ŝui k (f )| 2 Ŝuu k (f ) Ŝii k (f ) . (10)
The magnitude squared spectral coherence is then locally estimated around each block and averaged on an equivalent period of time controlled by the parameter α and Eq. ( 8).

In the two next sections, this algorithm is used and evaluated on simulated and experimental data.

Experimental approach

In this part, the time varying method of impedance estimation is tested on batteries. This study is performed on commercial lithium nickel manganese cobalt

Battery model and simulation

The goal of this section is to validate the previous approach thanks to a battery simulator. This simulator models the behavior of the battery at 25 • C.

To model the system behavior, an electric model with a constant phase element (CPE) [START_REF] Barsoukov | Impedance Spectroscopy -Theory, Experiment, and Applications[END_REF][25] is used. The impedance of such an element is given by:

Z cpe (w) = 1 Q(jw) p , ( 11
)
where Q is a real constant, j is a complex number such that j 2 = -1 and 0 ≤ p ≤ 1.

From Eq. ( 11), it is clear that if p = 1 the CPE behaves as a simple capacitor whereas if p = 0 it is considered as a resistor. Z thq (w The value of each parameter is obtained with a Levenberg-Marquardt [23][24] optimization algorithm applied to impedances estimated through narrowband FFT-EISs. These impedances are measured on a real battery under a constant discharge current of -0.

) = R 1 + jwL 1 + R 2 1 + jwR 2 C 1 + R 3 1 + R 3 Q(jw) p (12)

Simulation protocol

A variation of SoC from 90 % to 40 % is simulated by discharging the battery 325 with a constant polarization current of -0.5A which corresponds to a discharge ratio of C 4.4 . A PRBS current is continuously added to this polarization current in order to perform impedance estimation and tracking. The amplitude of this additive PRBS signal is set to 250 mA, which is sufficiently small for the battery to respect the linear assumption. In order to be in the same conditions as for the experimental results, the frequency band where the battery impedance is estimated is limited to 20 Hz -90 Hz. This band is chosen because of instrumental constraints which are detailed in the experimental protocol. This frequency bandwidth corresponds to a PRBS block of length 0.25 s. This block is then filtered thanks to a numeric low-pass filter with a cutoff frequency of 120 Hz to avoid abrupt variations in this additive current. The filtered block is then repeated continuously 31 680 times, i.e. during 2.2 h. The current and the voltage of the battery are, as in the experimental test, synchronously sampled at a sampling rate of 2 500 Hz. The simulated temperature is assumed constant at 25 • C during the discharge.

Experimental protocol

The battery is placed in an enclosure with controlled temperature of 25 • C. Then we consecutively apply the same block of filtered PRBS block, which was described in the protocol of the simulation, centered again around -0.5 A in order to apply global (Eq.(3b)) and exponential averaging and compare the obtained results. This low polarization current allows to minimize the SoC variation during the measurement time.

This filtered block is repeated consecutively 36 times, leading to a measurement process with a total duration of 9 s. During this process, the current and the voltage of the battery are synchronously sampled at a rate of 2500 Hz. The bandwidth studied begins at Fmin=20 Hz and closes at Fmax=90 Hz, due to instrumental constraints.

The chosen frequency band is also limited by the context of embedded applications. This clearly limits the measurement time and consequently the lower frequency at which the impedance can be correctly estimated. Indeed both the range of the sampling frequency and the number of points of an acquisition are limited. The lower Fmin is, the longer the SBPA block is. By choosing Fmax, the bit frequency Fbit of the PRBS block is defined as F bit = F max 0.4 [START_REF] Pintelon | System Identification -a Frequency Domain Approach[END_REF] and the sampling frequency Fe as at least F e = 10F bit. To raise Fmax, the sampling frequency should also be increased, yet for our experiment the sampling frequency is limited at most to 5 kHz. Moreover the higher the sampling frequency is, the larger the memory is.

The chosen frequency range is adapted to the study of the considered battery by an embedded measurement system. This frequency band corresponds approxi-

matively

Results and discussion

Simulation results

The time-varying impedance spectroscopy method developed in section 2 and summarized in Fig. 1 is applied to the simulated signals for the frequency band of interest (20 Hz -90 Hz). A Hann window and a forgetting factor α = 0.9 are chosen to estimate the necessary PSDs and CPSD with Eq. ( 7), the time-varying impedance with Eq. ( 9) and the time-varying magnitude squared spectral coherence with Eq. ( 10). This forgetting factor value corresponds to a local averaging over N = 19 blocks (see Eq. ( 8)) or equivalently a local averaging time of 4.75 s, and a response time at 80 % of 3.56 s (see Table . 1).

After this transient period, the algorithm is in steady state and the estimations can be considered as valid. On the frequency band of interest, the magnitude squared spectral coherence stays above 0.99 all along the simulation. This confirms that the battery can be considered as a LTI system during the equivalent averaging time, and that its electrochemical impedance can be correctly estimated. Fig. 4 presents the Bode diagram of the estimated impedance obtained with these settings for different SoC values during the whole discharge. After each decrease of 10 % of SoC, the estimated impedance is compared to its theoretical value computed thanks to Eq. ( 12). Fig. 4 shows that whatever the SoC value, the estimated impedance obtained 18 with the proposed method is very close to its theoretical value. The relative bias and standard deviation of the impedance modulus and phase estimation error averaged over the frequency band of interest are calculated for each SoC value.

During the whole discharge, these quantities stay below 1 % showing that the proposed method reaches good estimation and tracking performance. Moreover, a movie is available in supplementary information to further illustrate the tracking performance of the algorithm, and shows the real-time evolution of the estimated impedance during the whole discharge.

Finally, these simulation results show the good performance of the proposed method able to efficiently estimate and follow in real-time the temporal variations of a battery electrochemical impedance. In the next section, the proposed method is validated experimentally on a commercial battery.

Experimental results

During the whole measurement process, the SoC evolution is lower than 1% leading to the assumption that the battery can be considered as linear and time-invariant. This important assumption can be confirmed thanks to the magnitude squared spectral coherence defined by Eq. [START_REF] Blanke | Impedance measurements on leadacid batteries for state-of-charge, state-of-health and cranking capability prognosis in electric and hybrid electric vehicles[END_REF], where CPSD and PSD are estimated with α = 0.9 and an Hann window. Each time it stays higher than 0.99 in the frequency band of interest like in the simulation results. This high value confirms that during one estimation process, the battery behaves as a LTI system and its electrochemical impedance can be theoretically defined by Eq. ( 1) and estimated with Eq. ( 9). Fig. 5 shows the Nyquist plot results obtained, like the simulation results, with a Hann window and a forgetting factor α = 0.9 (corresponding to N=19 blocks following Eq. ( 8)) which induces a response time at 80 % of 3.56 s. We chose this value of α to converge before 9 s and to provide several impedance estimations after the convergence time. We can notice on the Nyquist plots that the global averaging estimation is very close to the estimation obtained through the FFT-EIS. However the exponential averaging method seems less accurate. This can be explained by the time window over which each estimator is averaged. The global averaging estimation is averaged over the whole data set (N = 36 blocks) leading to good estimation performance. On the contrary, the exponential averaging estimation is averaged over a smaller local window (α = 0.9, N = 19 blocks) whose size is set by using the forgetting factor α (the closer to 1 α is, the longer the time window is). The error curves shown in Fig. 5 represent the time evolution of the estimation performance indicator defined in the exponential (red curve) and global (green curve) cases. This indicator corresponds to the normalized root-meansquare error (RMSE Eq.( 13)) between a known reference value of the impedance Z ref (f ) and the k th estimated value of this impedance Z k (f ) described in Eq. [START_REF] Shin | Fundamentals of Signal Processing for Sound and Vibration Engineers[END_REF]. As expected, the error for the global averaging algorithm is always lower than the one obtained in the exponential averaging case. Moreover two steps are clearly visible on the red error curve: a first decaying part corresponding to the convergence time at 80 % of the algorithm (duration = 3.56 s for α = 0.9), and a second part after convergence where the estimate fluctuates around the correct impedance. We can see the same two steps on the film of the discharge on the simulator, after the convergence time, the estimation is really more accurate. These results emphasis the existence of a trade-off between the obtained convergence time (related to the tracking capabilities of the algorithm) and the final estimation error (related to the estimation performance of the algorithm). This trade-off can be managed through the value of the forgetting factor α. Indeed, a small forgetting factor leads to an algorithm able to follow strong variations of the estimated impedance over time, but with a large estimation error. On the contrary, a forgetting factor close to 1 leads to an algorithm able to precisely estimate the battery impedance, but unable to follow large time variations.

RMSE k = 100 f ∈B |Z ref (f ) -Ẑk (f )| 2 f ∈B |Z ref (f )| 2 (13)
To check the repeatability of the time-varying method, the experiment was repeated during 5 charges discharges cycles. The repeatability of the exponential method was evaluated and compared to the FFT-EIS performances.

In a first part, the aim is to quantify the repeatability of the FFT-EIS. For each SoC, 5 FFT-EISs estimations are averaged according to Eq. ( 14) and the standard deviation defined by Eq. ( 15) is calculated. Fig. 8 presents on a Bode diagram the averaged module and phase for several SoC.

Z eis m (f, SoC) = 1 5 i dech =5 i dech =1 Z eis i dech (f, SoC). (14)
Z eis std (f, SoC) = 1 4 i dech =5 i dech =1 |Z eis i dech (f, SoC) -Z eis m (f, SoC)| 2 . (15)
where Z eis i dech (f, SoC) is the impedance estimator for one SoC; Z eis m (f, SoC) corresponds to averaged impedance and Z eis std (f, SoC) its standard deviation. On the module part (on Fig. 8), it is very easy to distinguish each SoC, whereas it is more difficult on phase because of the inflection point around 25 Hz. The relative standard deviation average on the frequency bandwidth of the module is lower than the one of the phase, however in both cases these quantities are under 1.2 %, which show the good repeatability of this method. After having quantify the repeatability of the FFT-EIS, the tracking method is tested in a second part.

For each SoC, once convergence time is past, 5 estimations obtained with the exponential method are averaged by Eq. ( 16) and the standard deviation is calculated by Eq. [START_REF] Huet | A review of impedance measurements for determination of the state-of-charge or state-of-health of secondary batteries[END_REF].

Z m (f, SoC) = 1 5 i dech =5 i dech =1 Zi dech (f, SoC). ( 16
) Z std (f, SoC) = 1 4 i dech =5 i dech =1 | Zi dech (f, SoC) -Z m (f, SoC)| 2 . ( 17
)
where Zi dech (f, SoC) is the temporal average of impedance estimator for all the estimations after the convergence time for one SoC; Z m (f, SoC) corresponds to averaged impedance on 5 cycles and Z std (f, SoC) its standard deviation.

The results are presented on Bode diagram in Fig. 8 for several SoC. If on the module part it is again easy to follow the impedance evolution during the discharge, the phase presents more variations at some frequencies. Again an inflection point can be seen but with more difficulty than one in the FFT-EIS case because of the noise on the estimates.

The phase and module relative bias averages on the frequency bandwidth between the estimation by the proposed method and the FFT-EIS are calculated.

For each SoC the relative bias of the module (phase respectively) is under 1 % (respectively 2 %). The relative standard deviation average on the frequency bandwidth is under 0.1 % for the module and 1.1 % for the phase. These relative standard deviations are of the same order of magnitude than FFT-EIS.

The battery impedance was not estimated for extreme SoC values of 100 % and below 10 %. For a SoC of 100 %, the main reason is that measurement processes were carried on only during battery discharges. Under these conditions, the estimation of the battery impedance with a SoC of 100 % presents no practical interest and was not considered in the paper. For SoC values below 10 %, technical difficulties are encountered to estimate the battery impedance.

Indeed, the results obtained with the proposed method are compared to those obtained with FFT-EIS method. Unfortunately, no security system prevents the battery from deterioration at the end of the discharge during FFT-EIS. The impedance estimation for SoC 40 % to 10 % is supplied as supplementary information (Fig. B).

On this study over 5 discharges of the lithium-ion battery, the exponential method presents for each SoC low relative standard deviations and bias. The exponential method can so be considered as repeatable as the FTT-EIS method.

Conclusion

More and more applications use batteries. The development of BMS is needed to optimize their efficiency and dependability. With this in mind, estimate in real time the electrochemical impedance temporal evolution is a valuable tool. The method presented in this paper provides one solution to this open problem.

This method not only allows an accurate estimation of the battery impedance, but also a tracking of its temporal evolution. In other words, the developed method presents the major advantage to be evolutionary and to follow the temporal evolution of the impedance compared to classical estimation methods like FFT-EIS. It mainly relies on a simple exponential averaging strategy, and allows to manage the trade-off between its tracking and estimation performance through a single normalized parameter: the forgetting factor.

The chosen exponential averaging strategy can be recursively implemented with a low memory depth. Moreover, the use of a wideband signal such as PRBS can be generated with very simple electronics circuits. These two particularities clearly show that the proposed method is particularly relevant for embedded applications.

The proposed estimator was tested on both a simulated and real lithium ion battery, because this technology is widely used for electric or hybrid transports or connected objects. However, this method is not limited to lithium ion batteries and could also be applied to any other electrochemical systems such as primary battery, secondary battery (including aqueous technology or other) or fuel cells .

On the simulator, the exponential averaging method allows to track in real time the temporal variations of electrochemical impedance. The obtained results demonstrate the capabilities of the method to precisely track the battery impedance evolution over more than 50 % of SoC discharge.

On the real battery we checked the influence of the forgetting factor on the performance of estimation. A trade-off between the tracking and estimation per-

Figure 1 :

 1 Figure 1: Structure of the time-varying estimation of the battery electrochemical impedance.

Figure 2 :

 2 Figure 2: Equivalent electric model of the battery.

Fig. 2 and

 2 Fig. 2 and Eq. (12) present the equivalent electric model of the battery. Each element is attributed to various physicochemical phenomena in the battery and the measurement chain [25][26][27]. Thus, L 1 corresponds to the inductance of connexions, R 1 is the purely ohmic resistance of the battery and the connexions, R 2 is linked to the charge transfer resistance and C 1 to the double layer capacity, R 3 and Q, p model the diffusion of ions.

Fig. 3

 3 Fig. 3 clearly emphasis the variation of impedance during the discharge. It presents the Nyquist diagram of the modeled battery impedance for a SoC of 90 % and 40 %, from 0.1 to 1000 Hz in blue, obtained from Eq. (12) and Table. 2. The red curves correspond to the frequency band of interest (from 20 to 90 Hz) on which the time-varying method is tested in the experimental part.320

Figure 3 :

 3 Figure 3: Nyquist diagram of the modeled and measured impedance for a SoC of 90 % and 40 %.

  After a complete charge (25 • C, constant current of C 2 until 4.2V, constant voltage during 1h), we discharge the battery at a current of -0.5 A which corresponds to a rate of discharge of C 4.4 . The estimates obtained by the time-varying wideband frequency algorithm is compared to a known reference impedance called Z ref (f ). This reference is measured through an impedance spectroscopy thanks to a potentiotat VMP3 of Bio-Logic SAS with a booster 5A, associated with the EC-Lab software (signal amplitude of 200 mA, 1 Hz to 1 kHz, logarithm spacing of 10 measures per decade, 3 measures per frequency). Therefore, once 10 % of SoC is discharged, we first record a classical FFT-EIS in a galvanostatic mode at the polarization current. For instance, Fig. 3 presents the Nyquist diagram of this NMC battery obtained from the FFT-EIS for the SoC 90 % and 40 %.

  to the end of the first semicircle on the Nyquist diagram of the battery impedance. This area is closely related to the contributions of both ohmic and charge transfer resistance of the battery, which is an important parameter related to common state indicators such as the SoC and the SoH. Being able to estimate the battery impedance in this particular frequency band thus provides an efficient way to estimate this parameter and obtain information on the battery SoC and SoH. Indeed in the chosen frequency band, the battery impedance varies significantly with respect to the SoC, justifying the development of an evolutionary method. This point is highlighted by Fig. A in complementary information.

Figure 4 :

 4 Figure 4: Bode diagram of the estimated impedance during a discharge for different SoC values (top: modulus; bottom: phase).

Figure 5 :

 5 Figure 5: Estimated impedances between 20 Hz to 90 Hz by using the three methods (top:Nyquist diagram for exponential averaging α = 0.9 (red), global averaging (green), and the EIS method (blue), bottom: error curves (time evolution of the root-mean-square error)).

Figure 6 :

 6 Figure 6: Difference between the initial estimation and the first estimation after the convergence for the exponential case (left: 1 st estimation at time t=0.25 s, right: estimation after the convergence at time t=3.56 s, top: module, bottom: phase).

Fig. 6

 6 Fig. 6 underlines the difference of estimation accuracy between the first and after convergence temporal estimates of the battery impedance (both module and phase are represented). Once the convergent step of the exponential averaging algorithm is finished, the estimation performance becomes clearly more accurate, explaining the large decrease in the error curve of Fig. 5. The parameter which controls the length of the convergence time and the quantity of

Figure 7 :

 7 Figure 7: Influence of the forgetting factor (top: Nyquist diagram for α = 0.9 (red), α = 0.5(black), and the EIS method (blue), bottom: error curves for α = 0.5 (black), α = 0.9 (red), and the sliding window algorithm (green)).

Figure 8 :

 8 Figure 8: Bode diagram of the average on 5 discharges of impedance at 6 SoC (top: Module; bottom: Phase).

Table 1 :

 1 Main characteristics of the considered averaging strategies.

Table 2 :

 2 Values of the model parameters obtained for a SoC of 90 % and 40 %.

	5 A from a SoC of 90 % to 40 % every 10 % change. Dur-

ing this continuous discharge, parameters L 1 , C 1 , Q, p stay constant whereas

formance was determined. Then a repeatability study was carried out, showing that the proposed method finally appears as repeatable as the FFT-EIS method.

The proposed method provides not only an accurate and repeatable estimation of the impedance but also a tracking of its fluctuations over time. In further work, this algorithm will be embedded, applied on a test bench, and the developed method will be tested on batteries of different chemistries.