Giacomo Di

Mauro Mariani
email: mariani@mat.uniroma1.t

Giacomo Di Gesù
email: giacomo.di-gesu@enpc.fr

FULL METASTABLE ASYMPTOTIC OF THE FISHER INFORMATION

Keywords: 2010 Mathematics Subject Classification. 35P20, 28D20, 49J45, 60F10 Fisher Information, Metastability, Γ-convergence, Semiclassical spectral theory, Witten Laplacian

We establish an expansion by Γ-convergence of the Fisher information relative to the reference measure e -βV dx, where V is a generic multiwell potential and β → ∞. The expansion reveals a hierarchy of multiple scales reflecting the metastable behavior of the underlying overdamped Langevin dynamics: distinct scales emerge and become relevant depending on whether one considers probability measures concentrated on local minima of V , probability measures concentrated on critical points of V , or generic probability measures on R d . We thus fully describe the asymptotic behavior of minima of the Fisher information over regular sets of probabilities. The analysis mostly relies on spectral properties of diffusion operators and the related semiclassical Witten Laplacian and covers also the case of a compact smooth manifold as underlying space.

Introduction

The Fisher information of a probability measure µ relative to a reference measure m on a smooth manifold is given by the expression

I(µ|m) := 1 2 |∇ | 2 dm = 2 |∇ √ | 2 dm, (1.1)
assuming that := dµ dm exists and is sufficiently regular (see the precise Definition 3.1 below). This is a classical and ubiquitous object measuring the discrepancy between two measures and appearing in various fields as Statistics, Information Theory and Statistical Mechanics.

In the context of Statistical mechanics the reference measure usually appears in the form m β = e -βV dx, where V is a suitable potential describing the interaction in the model, β is a positive parameter, proportional to the inverse temperature and dx is the volume measure corresponding to the absence of interaction. The main interest often lies in the dynamical properties of the associated heat flow evolution featuring m β as stationary measure. Under suitable regularity assumptions on V it may be described by the Fokker-Planck equation

∂ t u = ∆u + β div(u∇V ) (1.2)
or, on a pathwise level, by the overdamped Langevin equation Ẋt = -∇V (X t ) + 2β -1 Ẇt ,

where W t is Brownian motion. The Fisher information may enter the description and analysis of a Statistical Mechanics model in various ways: through functional inequalities, prominently the Log-Sobolev and HWI inequality (see e.g. [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF] and references therein and in particular [START_REF] Menz | Poincaré and logarithmic Sobolev inequalities by decomposition of the energy landscape[END_REF] for the β → ∞ regime considered here); as rate functional in Large deviation principles describing fluctuations of the empirical occupation measure of (1.3) when t → ∞, [START_REF] Donsker | Asymptotic evaluation of certain Markov process expectations for large time (III)[END_REF]; or also as a tool to give an alternative construction of the dynamics (1.2) as gradient flow with respect to the relative entropy functional [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF] (see also the discussion at the end of this introduction).

1.1. The low temperature regime and metastability. In this paper we study the asymptotic behaviour of the Fisher information relative to m β = e -βV dx in the low temperature regime, that is β → ∞. When the potential V admits several local minima, in this regime energy traps are created which slowdown the dynamics (1.3) around these minima and produce metastability effects.

In the heuristic picture, there are several relevant time scales that feature a nontrivial dynamics: on short time scales, i.e. when t = t(β) grows slowly as a function of β, the system essentially follows the deterministic gradient flow with respect to the potential V . In particular, all critical points of V are absorbing for the limiting dynamics. On a longer time scale (see below), the process just lives on critical points; however on this scale stable critical points (that is, local minima of V ) are absorbing for the dynamics (once visited, the system sticks there).

When observed on even larger time scales, typically exponentially large in β [9, Chap. [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF].], even among the local minima distinctions become observable: on the long run, rare but sufficiently large fluctuations may occur that allow the process to climb the mountain pass which leads to another local minimum. This tunneling or metastable transition among local minima favors the passage from a minimum point to a deeper, energetically more convenient minimum point. Thus, on the long run, the deeper the local minimum, the more time the process will spend around it.

The reversible diffusion model (1.3) is a paradigmatic model for metastable phenomena, that show up in the dynamical behavior of a large variety of complex real world systems, and can be regarded as a mark of dynamical phase transitions. It is often a detailed theoretical and numerical analysis of metastability that allows for a proper upscaling of microscopic models -defined at the atomic or molecular level -to macroscopic descriptions catching experimental observations in materials science, biology and chemistry. We refer to Kramers' influential [START_REF] Kramers | Brownian motion in a field of force and the diffusion model of chemical reactions[END_REF], to [START_REF] Hänggi | Reaction-rate theory. Fifty Years after Kramers[END_REF] and to [START_REF] Berglund | Kramers' law: Validity, derivations and generalisations[END_REF][START_REF] Lelièvre | Two mathematical tools to analyze metastable stochastic processes[END_REF][START_REF] Gesù | Boris Nectoux Jump Markov models and transition state theory: the Quasi-Stationary Distribution approach[END_REF] for recent brief reviews on some of the mathematical techniques used in the analysis of metastability.

One expects the metastable behavior to be encoded somehow in the asymptotic behavior of the Fisher information. We discuss below why Γ-convergence is the natural tool to characterize this asymptotic behavior, in particular in the context of metastability and convergence of multi-scale reversible dynamics. However, for the sake of readability, we first shortly describe our main result in Section 1.2 below, then informally explain how Γ-convergence can be used as a solid theoretical framework for metastability in Section 1.3, in particular in the context of large deviations of the occupation measure. We quickly review the applications of Γ-convergence in multi-scale gradient flows in Section 1.4 (see also [START_REF] Peletier | From diffusion to reaction via Γ-convergence[END_REF] where Γ-convergence is exploited to upscale an essentially one-dimensional metastable Fokker-Planck equation). We recall the definition of Γ-convergence in Section 2.2 below, and refer to [START_REF] Braides | Γ-convergence for beginners[END_REF] for a systematic treatment. 1.2. Informal description of the main result. In this paper a full expansion in the sense of Γ-convergence of the Fisher information is established, which indeed features the complete cascade of metastable scales as informally described above. Using standard techniques, it is not hard to see that, under mild assumptions on the potential V , the rescaled Fisher information

I β := 1 β 2 I(•|e -βV dx) (1.4)
Γ-converges (in the weak topology of probabilities) as β → ∞ to the linear functional

I(µ) := 1 2 dµ(x)|∇V | 2 (x).
On the dynamical level this corresponds to the convergence of the diffusion (1.2)-(1.3) to the deterministic transport along the gradient of V . Note that Z -1 β e -βV dx is the unique minimizer of I β , provided Z β = exp(-βV ) < ∞. On the other hand, if V admits several critical points, I(µ) vanishes on any probability that is concentrated on such points.

This suggests that one may obtain a non-trivial Γ-limit also when multiplying I β by a suitable diverging constant, so that, according to the heuristic picture described above, the new limiting functional vanishes only on probability measures concentrated on local minima. This is indeed the case, and βI β also Γ-converges to a functional J with the mentioned properties (see Theorem 2.3 for its explicit characterization).

To capture the exponential scales associated with tunneling between local minima we iterate the same procedure: assuming that there is a unique global minimum point x 0 of V and labelling by {x 1 , . . . , x n } the other local minimum points of V , we investigate for each k = 1, . . . , n the Γ-limit of βe βW k I β , where W k is the mountain pass barrier which separates x k from a deeper minimum (see Section 2 for precise defintions). Under generic assumptions, we prove that for each k = 1, . . . , n also βe βW k I β Γ-converges to some J k , which can again be explicitly characterized (see Theorem 2.3) and involves the prefactors already appearing in the famous Kramers formula for metastable critical times [START_REF] Kramers | Brownian motion in a field of force and the diffusion model of chemical reactions[END_REF] (see also the Γ-convergence result in [START_REF] Peletier | From diffusion to reaction via Γ-convergence[END_REF]). Since J k (µ) = 0 for all k's iff µ = δ x0 , this amounts to a full expansion of I β by Γ-convergence (see [START_REF] Braides | Γ-convergence for beginners[END_REF]Chap. 1.10] and Section 2 of this paper). In other terms, as far as infima over closed and open sets are concerned (that is, in the sense of Γ-convergence)

I β (µ) ∼ I(µ) + 1 β J(µ) + n k=1 1 β e -β W k J k (µ). (1.5)
As a consequence of (1.5) and the Donsker-Varadhan Large Deviations Principle [START_REF] Donsker | Asymptotic evaluation of certain Markov process expectations for large time (III)[END_REF] we infer that, in the sense of large deviations, it holds (see Corollary 2.5)

P 1 T T 0 dsf (X s ) ∼ dµ(x)f (x), ∀f ∼ e -βT I(µ)+ 1 β J(µ)+ n k=1 e -β W k β J k (µ) . (1.6)
That is, we characterize the sharp asymptotic of the l.h.s. for every µ ∈ P(R d ).

To prove our main result (1.5) we take a spectral point of view and exploit results which were mainly developed in the context of the semiclassical spectral analysis of Schrödinger operators. Indeed, as highlighted by the expression on the right hand side of (1.1), the Fisher information can be seen as a Dirichlet form on L 2 (dm) and thus be studied from a spectral point of view by considering the corresponding selfadjoint operator. When m β = e -βV dx the latter is given by the diffusion operator

L β = ∆ -β∇V • ∇. (1.7)
It is well-known that a unitary transformation, sometimes called ground state transformation, turns the generator into a Schrödinger operator acting now on the flat space L 2 (dx). More precisely,

-e -1 2 βV L β e 1 2 βV = -∆ + β 2 4 |∇V | 2 - β 2 ∆V. (1.8)
The latter operator is a specific Schrödinger operator. This operator coincides with the restriction on the level of functions of the Witten Laplacian, an operator acting on the full algebra of differential forms, considered by Witten in his celebrated paper [START_REF] Witten | Supersymmetry and Morse Theory[END_REF]. Note that the small noise limit β → 0 turns now into a problem in semiclassical analysis of Schrödinger operators. In this paper we use two important results, established in this framework: the approximation of the low-lying spectrum of Schrödinger operators via harmonic oscillators sitting at the bottom of the wells [START_REF] Simon | Semi-classical analysis of low lying eigenvalues. Non-degenerate minima: asymptotic expansion[END_REF], and the fine analysis of the splitting of the exponentially small eigenvalues provided in [START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach[END_REF]. For simplicity we restrict here to the case of R d or of a compact Riemannian manifold as state space, but our arguments could be adapted without much difficulty to the case of a bounded domain with reflecting (Neumann) boundary conditions, as considered e.g. in [START_REF] Peletier | From diffusion to reaction via Γ-convergence[END_REF], by using [START_REF] Peutrec | Small eigenvalues of the Neumann realization of the semiclassical Witten Laplacian[END_REF] instead of [START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach[END_REF].

1.3. Metastable dynamical systems. Here we informally describe some motivations to approach metastability via Γ-convergence of the Fisher information. By the Birkhoff ergodic theorem, time-averages of an observable of an ergodic (deterministic or random) dynamical system converge in the long time limit to the average of the same observable with respect to the invariant measure m of the system. Assume that we have a family of ergodic systems indexed by some parameter β, that for the sake of simplicity we may take as real number. Under very weak and general conditions, see e.g. [START_REF] Kifer | Large deviations in dynamical systems and stochastic processes[END_REF], dynamical systems also satisfy a large deviation principle, that we informally write for each fixed β as

P β 1 T T 0 dsf (X s ) ∼ dµ(x)f (x), ∀f exp (-T E β (µ)) , T 1, (1.9)
where the rate E β (µ) ≥ 0 is a positive real valued function defined on probability measures 1 . For each fixed β, we assume that the system has good ergodic properties, i.e. E β (µ) = 0 iff µ = m β is the invariant measure 2 , so that (1.9) can be interpreted as a quantitative version of the Birkhoff theorem. Suppose now that, informally speaking, a limiting behavior takes place as we move the parameter β to some limit, say β → ∞. Then the invariant measure is 1 Notice that (1.9) makes sense both for random and deterministic systems, as in the latter case the initial distribution induces a probability measure on the trajectories of the system. In other words, for the deterministic systems, the l.h.s. of (1.9) should be interpreted as the measure of initial data such that T 0 f (Xs)ds ∼ dµ(x)f (x). 2 For the sake of simplicity we are assuming that at time 0 we start with distribution m β . This assumption is usually not needed in the random case, but crucial for deterministic systems.

expected to converge to some limiting measure, m β → m. Similarly one expects concentration of time averages to take place in the parameter β as well, namely

P β 1 T T 0 dsf (X s ) ∼ dµ(x)f (x), ∀f exp (-T a β I(µ)) , T 1, β 1 (1.10)
where a β is a sequence of real numbers converging to +∞. It is a general fact, see e.g. [START_REF] Mariani | A Γ-convergence approach to large deviations[END_REF], that the functional I in (1.10) is nothing but the Γ-limit of 1 a β E β . If ergodic properties of the system hold uniformly in β, one would expect that

I(µ) = 0 iff µ = m = lim β m β .
On the other hand, we say that the system exhibits a metastable behavior (on the scale (a β )) in the limit β → ∞, if the function I(µ) = Γ -lim β 1 a β E β , vanishes on a set strictly larger than {m}. In other words, if there is a 'small-but-not-so-small' probability that, for T large and β large, the time average of observables does not behave as the spacial average w.r.t. the limiting invariant measure. In such a case, (1.10) states in particular that one needs to consider time scales much longer than a β to observe convergence of the system to the invariant measure. If the system features a metastable behavior, (1.10) is not completely satisfactory. Indeed, one knows a priori that the l.h.s. should vanish if µ = m (as we sent T → ∞ before β, and assumed m β → m), but the r.h.s. does not catch the sharp asymptotic if I(µ) = 0 (and µ = m). In other words, if µ = m and I(µ) = 0 (1.10) just states that the l.h.s. vanishes on a slower scale than exp(-a β T ) but does not quantify this scale.

In this case, one would expect a further non-trivial expansion in the r.h.s. of (1.10), say for T 1, β 1

P β 1 T T 0 dsf (X s ) ∼ dµ(x)f (x), ∀f exp -T a β I(µ) + 1 b β J(µ) , (1.11)
where b β → ∞ and J(µ) is such that J(m) = 0, J(µ) = ∞ if I(µ) > 0, and thus J(µ) gives the asymptotic value of the l.h.s. when I(µ) = 0 and µ = m. Actually, it may happen that for such µ's, the l.h.s. in (1.11) vanishes at different rates, so that a further expansion of the exponent in the r.h.s. is necessary, until the exponential behavior of the l.h.s. at each and every point in {I(µ) = 0} is characterized, and a complete quantitative version of the ergodic convergence of averaged observables is recovered.

The discussion informally carried over above can be made very precise, since both large deviations and Γ-convergence always hold along subsequences, and thus (1.9) can actually be used as a rigorous definition of E β , and the existence of a nontrivial development by Γ-convergence [4, Chap. 1.10] (or equivalently of multiple large deviations principles3 as in (1.11)) as a general definition of metastability. However we refrain from giving a too abstract formulation of these statements that are not instrumental to our results below.

1.4. Multiscale gradient flows. The relative entropy H(µ) of the probability measure µ w.r.t. the reference measure m = e -βV dx is defined as

H β (µ) := dm log , µ = m (1.12)
It is a well-known fact [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF] that the Fokker-Planck equation

∂ t u = ∆u + βdiv(u∇V ) (1.13)
can be written as a gradient flow of the relative entropy H with respect to the 2-Wasserstein distance on probability measures on R d . Namely a curve (µ β (t)) of probability measures with bounded second moments on R d satisfies (1.13) iff

H β (µ β (0)) ≥ H β (µ β (t)) + 1 2 t 0 ∂ s µ β (s) 2 ds + 1 2 t 0 |∇ W2 H β | 2 (µ β (s))ds (1.14)
where • is the tangent norm on the 2-Wasserstein and, most importantly for us, the metric gradient |∇ W2 H β | 2 coincides with βI β . Following the approach in [START_REF] Sandier | Gamma-convergence of gradient flows with applications to Ginzburg-Landau[END_REF], and dividing both sides of (1.14) by β, one wants to pass to the limit β → ∞.

Then, using the Γ-limit result I β → I in Theorem 2.3, one easily gathers that any limit point (µ(t)) of (µ β (t)) satisfies

V dµ(0) ≥ V dµ(t) + 1 2 t 0 ∂ s µ(s) 2 ds + 1 2 t 0 |∇V | 2 dµ(s)ds, (1.15)
provided the initial data satisfies

1 β H β (µ β (0)) → V dµ(0).
In other words, one recovers the simple fact that the solution u β to (1.13) satisfies u β (t/β, x) → u(t, x), where u solves the dual transport equation

∂ t u = div(u∇V ).
(1.16)

Of course, if the inital data concentrates on critical points of V , (1.16) becomes trivial, namely ∂ t u = 0. Then the development by Γ-convergence gives the sharp asymptotic of u β (t/β 2 ) exactly in these cases. It turns out that, when the initial condition concentrates on local minima, even this limit is trivial, and still by the Γ-development one can establish the limit of u β (t/(βe βW k ), x) where the W k are appropriate constants defined in A.5 below. We do not pursue these problems here, since we plan a detailed study in the non-reversible case (namely when ∇V in (1.16) is replaced by a generic vector field) in a follow-up paper. See [START_REF] Peletier | From diffusion to reaction via Γ-convergence[END_REF] for a detailed discussion on a similar model.

1.5. Plan of the paper. In Section 2 we introduce some definitions and state the main result. In Section 3, we shortly review some basic properties of the Fisher information functional. In Section 4 we cosnider suitable quasimodes associated to the generator L β . In Section 5 we prove the Γ-convergence and equicoercivity results for the Fisher information. In Section 6 and Section 7 we prove the development by Γ-convergence of the Fisher information respectively under inverse power and exponential rescaling.

Main result

2.1. Basic notation. Hereafter (M, g) is a smooth d-dimensional Riemannian manifold without boundary and with metric tensor g. We assume that either M is compact or that M = R d and g is the canonical Euclidean tensor. Ω is the set of smooth, compactly supported 1-forms over M . For ψ ∈ Ω, ψ 2 is understood as ψ, ψ g , while ∇, ∇•, Hess and ∆ are the covariant gradient, the divergence, the Hessian and the (negative) Laplace-Beltrami operators on M . P(M ) is the set of Borel probability measures on M , and it is naturally equipped with the weak topology of probability measures. If M is compact this is the weakest topology such that the map µ → dµ f is continuous for all f ∈ C(M ). If M is not compact, the definition is slightly more involved, as one should require f to be bounded and uniformly continuous with respect to any totally bounded distance on M (this topology is independent on such a distance), see [START_REF] Stroock | Probability Theory. An analytic view[END_REF]Chap. 3.1]. In any case, P(M ) is a Polish space (that is a completely metrizable, separable topological space), and it's compact if M is compact.

2.2. Γ-convergence. In this section we briefly recall the basic definitions related to Γ-convergence, see [START_REF] Braides | Γ-convergence for beginners[END_REF]Chap. 1]. For X a Polish space (we will always consider X = P(M )) and (H β ) a family of lower semicontinuous functions on X indexed by the directed parameter β (we will always consider β > 0), one defines

(Γ -lim β H β )(x) := inf lim β H β (x β ), {x β } ⊂ X : x β → x , Γ -lim β H β (x) := inf lim β H β (x β ), {x β } ⊂ X : x β → x .
(2.1)

Whenever Γ -lim β H β = Γ -lim β H β =: H we say that H β Γ-converges to H in X.
Equivalently, H β Γ-converges to H iff for each x ∈ X:

-for any sequence x β → x it holds lim β H β (x β ) ≥ H(x); -there exists a sequence

x β → x such that lim β H β (x β ) ≤ H(x). (H β ) is equicoercive if for each M ∈ R the set lim β H β ≤ M is precompact in X.
If H β is equicoercive and Γ-converges to H, then A := argmin(H) contains each limit point of argmin(H β ). Denote by A 0 the set of such limit points. If A \ A 0 is nonempty, there exists a diverging sequence (a

β ) such that the functional

H (1) β := a (1)
β (H β -inf x∈X H β (x)) (2.2)
admits a non-trivial (that it having somewhere a value different from 0 and ∞) Γ-liminf. Notice that H

(1) β inherits from H β the property of being equicoercive, and that (Γ -lim β H

(1)

β )(x) = (Γ -lim β H (1) β )(x) = 0 if x ∈ A 0 and (Γ -lim β H (1) β )(x) = (Γ -lim β H (1) β )(x) = +∞ if x ∈ A. If H (1
) β admits a Γ-limit H (1) then we say that the development by Γ-convergence

H β ∼ H + 1 a (1) β
H (1) (2.3) holds. Iterating the procedure with sequences such that a

(k) β /a (k+1) β
→ 0, we say that the development by Γ-convergence is full if it holds

H β ∼ H + k 1 a (k) β H (k) (2.4)
and for each x ∈ A \ A 0 there exists k such that

H (k) (x) ∈ (0, ∞).
Equicoercivity and Γ-convergence of a sequence (H β ) imply an upper bound of infima over open sets, and a lower bound of infima over closed sets, see e.g. [START_REF] Braides | Γ-convergence for beginners[END_REF]Prop. 1.18], and it is a relevant notion of variational convergence for the problems discussed in the introduction. A full development by Γ-convergence then gives a sharp asymptotic of each of such infima.

2.3.

Assumptions on the potential. The potential V is a given real-valued function and we consider the following assumptions.

A.1 V ∈ C ∞ (M ) is a Morse function. Namely the Hessian of V is non- degenerate at any critical point of V . A.2 In the case M = R d (M is compact otherwise), it holds for some β 0 ≥ 0 lim x→∞ inf β≥β0 ( 1 4 |∇V | 2 - 1 2β ∆V )(x) = +∞.
(2.5)

A. [START_REF] Bovier | Metastability in reversible diffusion processes II. Precise asymptotics for small eigenvalues[END_REF] In the case M = R d (M is compact otherwise), it holds for some β 0 ≥ 0 and all β ≥ β 0

Z β := e -βV (x) dx < +∞. (2.6) Remark 2.1. If M = R d , A.1-A.3 imply that lim |x|→+∞ V (x) = +∞ (see [11, Proposition 2.2]). Thus, in both cases of M compact and M = R d , the set {V ≤ c} is compact for all c ≥ inf x∈M V (x).
Remark 2.2. Assumption A.2 says that, after the ground state transformation of the diffusion generator L β and division by β 2 (see (1.7),(1.8)), the Schrödinger potential grows to infinity at infinity, uniformly in β.

Let ℘ be the set of critical points of V , namely z ∈ ℘ iff ∇V (z) = 0. By A.1-A.2, ℘ is finite. Let ℘ 0 ⊂ ℘ be the set of local minima of V and define W : ℘ 0 → (0, ∞] as

W (x) := inf y =x : V (y)≤V (x) inf γ∈Γ(x,y) sup t∈[0,1] V (γ(t)) -V (x), (2.7)
where Γ(x, y) is the set of continuous curves connecting x and y in time 1. In other words W (x) is the lowest mountain pass to climb when going from x to a deeper minimum of V .

A. [START_REF] Braides | Γ-convergence for beginners[END_REF] For each x ∈ ℘ 0 such that W (x) < +∞ there exists a unique point x ∈ ℘ such that the two following conditions hold:

(i) V (x) = V (x) + W (x).
(ii) x and x lie in the same connected component of the compact set {y ∈ M : V (y) ≤ V (x) + W (x)}. In other words, there is a unique saddle point x such that all the optimal curves γ in the variational problem (2.7) pass through x. A.5 W (x) = W (y) whenever x = y, x, y ∈ ℘ 0 .

By A.5 we can label ℘ 0 = {x 0 , . . . , x n } by requiring

W n < . . . < W 1 < W 0 = +∞, (2.8)
where we use the shorthand notation W i := W (x i ). Note that x 0 is the unique global minimizer of V and xi exists for i = 0 by A.4.

2.4.

The Γ-development theorem. For β ≥ β 0 define the reference measure m β on M as

dm β (x) = e -βV (x) dx, (2.9)
where dx is the Riemannian volume on M . The Fisher information at inverse temperature β > 0 is the functional (see Definition 3.1 and (3.6) for further details)

I β : P(M ) → [0, +∞], I β (µ) := 1 2β 2 dm β |∇ | 2 if dµ dm β = ; +∞ otherwise.
(2.10)

We are interested in the variational convergence of I β in the low temperature regime, namely when β → ∞. In order to describe our main result we need to define further notation.

For z ∈ ℘, we denote by (ξ i (z)) i=1,...,d the eigenvalues of Hess V (z), labelled by ordering ξ 1 (z) ≤ ξ 2 (z) ≤ . . . ≤ ξ d (z). Then define

ζ : M → [0, +∞] ζ(z) := d i=1 |ξ i (z)| -ξ i (z) = 2 i : ξi(z)<0 |ξ i (z)| if z ∈ ℘; +∞ otherwise.
(2.11)

For k ≥ 1 and x k ∈ ℘ 0 , the saddle point xk defined in A.4 is easily seen to satisfy ξ 1 (x k ) < 0 and ξ i (x k ) > 0 for i ≥ 2. Recalling the labelling (2.8) we define for k = 1, . . . , n

η k : M → [0, +∞] η k (x) :=          0 if x ∈ {x 0 , . . . , x k-1 }; |ξ1(x k )| π det(Hess V )(x k ) | det(Hess V )(x k )| = |ξ1(x k )| d i=1 ξi(x k ) π 2 d i=2 ξi(x k ) if x = x k ; +∞ otherwise.
(2.12)

The functions ζ and η k are lower semicontinuous and coercive. Their expression naturally appears when studying the asymptotic behaviour of eigenvalues of L β (see Theorem 4.1 and 4.5 for this interpretation). Our main result is the following.

Theorem 2.3. Assume A.1-A.5. Then (I β ) is an equicoercive sequence of lower semicontinuous functionals. Moreover Γ -lim

β I β = I, (2.13)
Γ -lim

β βI β = J, (2.14)
Γ -lim

β (β e βW k I β ) = J k for k = 1, . . . , n, (2.15)
where the lower semicontinuous, coercive functionals I, J, J k : P(M ) → [0, +∞] are defined as

I(µ) := dµ(x) |∇V | 2 (x).
(2.16)

J(µ) = dµ(x) ζ(x) = y∈℘ α y ζ(y) if µ = y∈℘ α y δ y ; +∞ otherwise.
(2.17)

J k (µ) = dµ(x) η k (x) = α k η k (x k ) if µ = k j=0 α j δ xj ; +∞ otherwise. (2.18) Remark 2.4. It holds inf µ I(µ) = inf µ J(µ) = inf µ J k (µ) = 0.
The number of minimizers of the functionals in the sequence I, J, J n , . . . , J 1 decreases from |℘| to 1. More precisely we have

• I(µ) = 0 iff J(µ) < +∞, namely iff µ is concentrated on ℘. • J(µ) = 0 iff J n (µ) < +∞, namely iff µ is concentrated on ℘ 0 . • for k = 1, . . . , n-1, J k+1 (µ) = 0 iff J k (µ)
< +∞, and finally J 1 (µ) = 0 iff µ = δ x0 (that is the unique limit point of m β , the minimizer of I β ). In other words, I β admits the following full development by Γ-convergence, see [4, Chap. 1.10], as β → ∞:

I β (µ) ∼ I(µ) + 1 β J(µ) + n k=1 1 β e -β W k J k (µ). ( 2

.19)

In particular the asymptotic behavior of infima of I β on closed and open subsets is characterized by a finite and non-zero leading order, see [4, Theorem 1.18].

Corollary 2.5. Consider the diffusion process X β defined through

Ẋβ t = -∇V (X β t ) + 2β -1 Ẇt X β (0) = x 0 . (2.20)
Equivalently, let X β be the Markov process with generator β -1 ∆ -∇V • ∇, and starting at x 0 ∈ M . For β, T > 0, consider the empirical measure

θ β, T := 1 T T 0 δ X β t dt ∈ P(M ). ( 2

.21)

In other terms, θ β,T is the random probability measure on M such that for all test functions ϕ it holds dθ β,T (ϕ) = 1 T T 0 ϕ(X β t )dt. Then, in the sense of large deviations [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF], and in the limit T >> β → +∞, θ β,T satisfies (independently of x 0 )

P(θ β,T ∼ µ) ∼ e -βT I(µ)+ 1 β J(µ)+ n k=1 e -β W k β J k (µ) . (2.22)
More precisely, for all I-regular (respectively J-regular and J k -regular) sets A ⊂ P(M ) it holds

lim β→∞ lim T →∞ 1 β T log P(θ β,T ∈ A) = -inf µ∈A I(µ), (2.23)
lim β→∞ lim T →∞ 1 T log P(θ β,T ∈ A) = -inf µ∈A J(µ), (2.24) lim β→∞ lim T →∞ e βW k T log P(θ β,T ∈ A) = -inf µ∈A J k (µ). ( 2

.25)

Proof. It is not hard to check that, under A.1-A.3, the hypotheses in [START_REF] Donsker | Asymptotic evaluation of certain Markov process expectations for large time (III)[END_REF] are satisfied, and thus for each fixed β > 0 the empirical measure θ β, T satisfies a large deviations principle with rate β I β as T → +∞. As well known and easy to prove, the rate function of the large deviations with speed a β T , for the directed familiy θ β, T as T → ∞ and β → ∞ is then given by the Γ-limit of 1 a β βI β , see [START_REF] Mariani | A Γ-convergence approach to large deviations[END_REF]. Thus the statement follows from Theorem 2.3 when using respectively a β = β, a β = 1 and a β = e -βW k .

The Fisher information

In this section we shortly collect some basic facts concerning the Fisher information. We only sketch the proofs here since we just restate some facts in what is a well understood framework. In this section we omit the dependence on β for the sake of readability. Definition 3.1. Let V ∈ C 0 (M ). The Fisher information I : P(M ) → [0, +∞] induced by V is defined as

I(µ) = sup ψ∈Ω dµ e V ∇ • ψ -1 2 e 2V ψ 2 . (3.1)
Define the measure m ∈ P(M ) as dm(x) = e -V (x) dx.

Remark 3.2. As a supremum of linear and continuous functionals, I is convex and lower semicontinuous. Moreover if V ∈ C 1 (M ), the change of variable ψ → -e -V ψ shows

I(µ) = sup ψ∈Ω dµ ∇V • ψ -∇ • ψ -1 2 |ψ| 2 = 1 2 dm |∇ | 2 if dµ = dm and ∇ log ∈ L 2 (µ) +∞ otherwise = dx 2|∇h| 2 + h 2 |∇V | 2 2 -∆V if dµ = h 2 dx and h ∈ H 1 (M ) +∞ otherwise. (3.2)
Ĩ(µ) = 1 2 dm |∇ | 2 if µ ∈ P(M ) +∞ otherwise. (3.4)
Proof. Since changing V to V + log Z does not change I, here we assume Z = 1. We need to prove that for each µ there exists µ n ∈ P(M ) such that µ n → µ and lim n Ĩ(µ n ) ≤ I(µ). Otherwise said, that P(M ) is I-dense in P(M ). Setting

µ n = (1 -1 n )µ + 1 n m (3.5)
one easily gathers µ n → µ and I(µ n ) ≤ I(µ). Therefore it's enough to prove that P(M ) is I-dense in the set of µ ∈ P(M ) with density bounded away from 0. This is immediate from the density of smooth functions in weighted Sobolev spaces.

In view of Remark 3.2, the functional I β , already introduced in (2.10) and appearing in our main theorem, can be equivalently defined as

I β (µ) := β -2 I βV , (3.6)
where I βV is the Fisher information induced by βV as defined in Definition 3.1.

Spectral analysis of the generator

In this section we denote by •, • the inner product in L 2 (m β ). Consider the operator

L β f := ∆f -β∇V • ∇f (4.1)
defined for f ∈ C ∞ c (M ). By A.2 and standard results, it uniquely extends to a self-adjoint operator in L 2 (m β ) with compact resolvent. We still denote by L β such an extension.

S λ := (z, n) ∈ ℘ × N d : λ z,n = λ . (4.3)
The following theorem conveniently resumes in our setting some of the main results obtained in the paper [START_REF] Simon | Semi-classical analysis of low lying eigenvalues. Non-degenerate minima: asymptotic expansion[END_REF] in a more general framework of semiclassical Schrödinger operators. Throughout this Subsection 4.1 it is sufficient to assume just A.1 and A.2. Theorem 4.1 (Bounded eigenvalues). Fix Λ > 0. There exists ε Λ > 0 such that for each ε ∈ (0, ε Λ ) the following holds: there exists β ε > 0 so that for β ≥ β ε Spec(-

1 β L β ) ∩ [0, Λ] ⊂ λ : S λ =∅ (λ -ε, λ + ε) ∩ [0, Λ]. (4.4)
Moreover, for each λ ∈ [0, Λ] such that S λ = ∅, the cardinality of Spec(-1 β L β )∩(λε, λ+ε) equals, when each eigenvalue is counted with its multiplicity, the cardinality of S λ .

This theorem states that in the limit β → ∞, counting multiplicity, the spectrum of -1 β L β is well-aproximated by the spectrum of the direct union of suitable Harmonic oscillators. It is necessary to fix a treshold Λ, since this harmonic approximation is not uniform for Λ → ∞.

For λ ≥ 0 and ε > 0 denote by

P β,ε,λ := 1 (λ-ε,λ+ε) (-1 β L β ) (4.5)
the spectral projection of -1 β L β associated to the interval (λ -ε, λ + ε). Proposition 4.2. There exists an orthonormal base (Ψ β,z,n ) z∈℘,n∈N d of L 2 (m β ) such that for every ε > 0, z ∈ ℘ and n ∈ N d the following holds: there exists a β ε,z,n > 0 such that for β > β ε,z,n and for every ϕ

∈ C b (M ), Ψ β,z,n ∈ Range P β,ε,λz,n , (4.6) lim β Ψ β,z,n Ψ β,z ,n ϕ dm β = ϕ(z) if z = z and n = n 0 otherwise . (4.7)
We stress that eigenbases of L β in general fail to satisfy (4.6)-(4.7), due to possible resonances (that is λ z,n = λ z ,n for some z = z and n, n ). We prove the proposition after the preliminary Lemma 4.3.

Lemma 4.3. For each z ∈ ℘, n ∈ N d there exists Ψβ,z,n ∈ C 2 c (M ) such that for every ϕ ∈ C b (M ) lim β |(-1 β L β -λ z,n ) Ψβ,z,n (x)| 2 dm β = 0, (4.8) lim β Ψβ,z,n Ψβ,z ,n ϕ dm β = ϕ(z) if z = z and n = n 0 otherwise . (4.9)
Proof. First consider first the case M = R d . Define for β > 0 the differential operator (see also (1.8))

A β := -∆ + 1 4 β 2 |∇V | 2 - 1 2 β∆V (4.10) and let for z ∈ ℘ H β,z := -∆ + β 2 Q z -βC z , (4.11)
where

Q z (x) = 1 4 Hess 2 V (z)(x -z) • (x -z), C z = 1 2
Trace(Hess V (z)). (4.12)

The operator H β,z extends to a self-adjoint operator in L 2 (dx), which is nothing but a harmonic oscillator shifted by the constant βC z . In particular its spectrum is given by (βλ

√ Z β,z,n
, where

χ z ∈ C ∞ c (M ; [0, 1]
) is a smooth cutoff function such that χ ≡ 1 on a ball B r (z) centered at z of radius r, and χ ≡ 0 on B c 2r (z), where r > 0 is chosen sufficiently small (in particular such that the χ z 's have pairwise disjoint supports).

It follows from Taylor expansion of the Schrödinger potential

1 4 β 2 |∇V | 2 -1 2 β∆V that, as x → z, A β = H β,z + O β 2 |x -z| 3 + O (β|x -z|) . (4.15)
Thus one obtains

(A β -βλ z,n ) Θβ,z,n (x) 2 dx = 1 Z β,z,n Θ β,z,n ∆χ z + 2∇χ z • ∇Θ β,z,n + O β 2 |x -z| 3 + O (β|x -z|) χ z Θ β,z,n 2 dx.
(4.16)

Hence, using the explicit form of the Θ β,z,n 's as Hermite functions centered at z, it follows by Laplace asymptotics that and notice that

1 β 2 (A β -βλ z,n ) Θβ,z,n (x) 2 dx = O(β -1 ). ( 4
|(-1 β L β -λ z,n ) Ψβ,z,n (x)| 2 dm β = 1 β 2 |(A β -βλ z,n ) Θβ,z,n (x)| 2 dx
and Ψβ,z,n Ψβ,z ,n ϕ dm β = Θβ,z,n Θβ,z ,n ϕ dx.

Therefore (4.8) and (4.9) follow from (4.17) and (4.18). The proof in the case of M a compact manifold follows from a straightforward adaptation of the previous arguments: since the Laplace-Beltrami operator in coordinates takes the form i,j

g ij ∂ i ∂ j + ∂ j g ij + g ij ∂ j log det g ∂ j , (4.19)
using in particular local coordinates {y = y (z) } around each z ∈ ℘ such that the metric tensor is the identity in 0, leads (instead of (4.15)) to the local estimate

A β = H β,z + i [O (|y|) ∂ i ∂ i + O(1)∂ i ] + O β 2 |y| 3 + O (β|y|) . (4.20)
Chosing r sufficiently small such that the support of χ z is contained in the corresponding local chart {y (z) }, one can argue as before through Laplace asymptotics and show that the additional terms appearing in (4.20) give again an O(β) contribution.

Remark 4.4. Let T be a bounded operator on a Hilbert space H and let u, v ∈ H be normalized. Then

| v, T v -u, T u | ≤ 2 T 1 -u, v 2 + 1 -u , v 2 . (4.21) Proof. Write v = u, v u + (v -u, v u) and notice v -u, v u H = 1 -u, v 2 . Using the inequality √ 1 -a 2 √ 1 -a 2 + √ 1 -aa ≤ √ 1 -a 2 + √ 1 -a 2
for a = u, v and a = u , v one gets the result by straigthforward linear algebra.

Proof of Proposition 4.2. We use the same notation as in Theorem 4.1. Fix Λ > 0 large enough. For λ < Λ such that S λ = ∅, for (z, n) ∈ S λ take Ψβ,z,n as in Lemma 4.3. Define Ψβ,z,n := P β,εΛ,λ Ψβ,z,n as the spectral projection of Ψβ,z,n corresponding to the interval (λ

-ε Λ , λ + ε Λ ), see (4.5). Now for (z, n), (z , n) ∈ S λ , since P β,εΛ,λ is a projector, Ψβ,z,n , Ψβ,z ,n L 2 (m β ) = Ψβ,z,n , Ψβ,z ,n L 2 (m β ) -Ψβ,z,n -Ψβ,z,n , Ψβ,z ,n -Ψβ,z ,n L 2 (m β ) . (4.22)
The definition of Ψβ,z,n readily implies

Ψβ,z,n -Ψβ,z,n L 2 (m β ) = (I -P β,εΛ,λ ) Ψβ,z,n L 2 (m β ) ≤ 1 ε Λ (-1 β L β -λ z,n ) Ψβ,z,n L 2 (m β ) . (4.23)
In view of (4.8) this vanishes as β → ∞. Therefore the last term in (4.22) vanishes as β → ∞ and (4.9) now gives

lim β Ψβ,z,n , Ψβ,z ,n L 2 (m β ) = 1 if z = z and n = n 0 otherwise. ( 4
H = L 2 (m β ), v = Ψ β,z,n , v = Ψ β,z ,n , u = Ψβ,z,n , u = Ψβ,z ,n and T ≡ T ϕ the pointwise multiplication by ϕ to get Ψ β,z,n , ϕΨ β,z ,n L 2 (m β ) -Ψβ,z,n , ϕ Ψβ,z ,n L 2 (m β ) ≤ 2 ϕ ∞ 1 -Ψβ,z,n , Ψ β,z,n 2 L 2 (m β ) + 1 -Ψβ,z ,n , Ψ β,z ,n 2 L 2 (m β ) . (4.26)
Thus we obtain (4.7) from Lemma 4.3-(ii) and (4.25).

Thus the (Ψ β,z,n ) (z,n)∈S λ are an orthonormal base of Range(P β,εΛ,λ ) satisfying (4.7). Since this holds for any λ such that S λ is not empty and L 2 (m β ) is a direct sum of eigenspaces of L β , the (Ψ β,z,n ) (z,n) are an orthonormal base of L 2 (m β ) as λ runs over the positive real numbers such that S λ is nonempty. 4.2. Eigenfunctions of the split spectrum. We denote by β,0 , . . . , β,n the first n+1 eigenvalues of -1 β L β counted with multiplicity. It follows from the smoothness of V , A.2, A.3 and general principles that β,0 equals zero for every β and is simple with constant eigenfunction Φ β,0 ≡ 1/ Z β . Moreover, with the notation of the previous Subsection 4.1, λ x,0 = 0 whenever x ∈ ℘ 0 is a local minimum, and λ z,n > 0 otherwise. It follows then from Theorem 4.1 that there exist ε > 0 and

β 0 > 0 such that for each β ≥ β 0 Spec(-1 β L β ) ∩ [0, ε] = { β,0 = 0, β,1 , . . . , β,n } . (4.27)
It is not difficult to see that the eigenvalues β,1 , . . . , β,n are indeed exponentially small in β. A much stronger statement, giving the precise leading behaviour of the β,k 's in a general setting encompassing hypotheses A.1-A.5, was given in [START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach[END_REF], (see also [START_REF] Bovier | Metastability in reversible diffusion processes II. Precise asymptotics for small eigenvalues[END_REF] and [START_REF] Eckhoff | Precise asymptotics of small eigenvalues of reversible diffusions in the metastable regime[END_REF]). The following theorem conveniently resumes in our framework a weak version of the main result in [START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach[END_REF].

Theorem 4.5 (Low-lying eigenvalues). For k = 1, . . . , n the eigenvalue β,k is simple, admits a normalized eigenfunction Φ β,k ∈ C ∞ (M ) and satisfies, for η k as in (2.12), lim

β e βW k β,k = 1 2 η k (x k ). (4.28)
Lemma 4.6. For k = 1, . . . , n there exists

Φβ,k ∈ C ∞ c (M ) with Φβ,k L 2 (m β ) = 1, such that (i) the probability measure Φ2 β,k dm β converges to δ x k in P(M ). (ii) lim β→∞ Φβ,k , Φ β,k 2 = 1 .
Proof. In this proof the W k 's are defined as in A.4-A.5, and we set for convenience W n+1 = 0 and β,n+1 = ε, with ε as in (4.27) (notice that this is not an eigenvalue of -

1 β L β ). Fix δ > 0 such that δ < min k=1,...,n (W k -W k+1 ). (4.29) Let B k ⊂ M , respectively B k,δ , be the connected component of V -1 ((-∞, V (x k ))), respectively V -1 ((-∞, V (x k ) -δ)), containing x k .
The B k 's are precompact by Remark 2.1, and the closure of B k,δ is contained in B k . In particular there exists

χ k ∈ C ∞ c (M ; [0, 1]) such that χ k ≡ 1 on B k,δ and χ k ≡ 0 on M \ B k . For k = 1, . . . , n define Φβ,k = χ k / Z β,k where Z β,k = e -βV (x) χ 2 k (x) dx. ( 4
x : |∇χ k (x)|>0 V (x) ≥ V (x k ) -δ = V (x k ) + W k -δ. (4.34) Denote now by P β,k := 1 [ β,k+1 ,+∞) (-1 β L β ) the spectral projection of -1 β L β asso- ciated to the interval [ β,k+1 , +∞). By spectral decomposition P β,k Φ = Φ - k j=0 Φ, Φ β,j Φ β,j , ∀Φ ∈ L 2 (m β ). ( 4

.35)

The Markov inequality for spectral projections, an integration by parts and (4.34) yield

P β,k Φβ,k 2 L 2 (m β ) ≤ 1 β,k+1 Φβ,k , -L β Φβ,k = 1 β,k+1 Z β,k |∇χ k | 2 (x)e -βV (x) dx ≤ C k exp(-β(V (x k ) + W k -δ)) β,k+1 Z β,k , (4.36)
where Taking the scalar product with Φβ,k in both sides of (4.35) calculated for Φ = Φβ,k , gives

C k = |∇χ k | 2 (x)dx < ∞.
1 -Φβ,k , Φ β,k 2 = Φβ,k , P β,k Φβ,k + k-1 j=0 Φβ,k , Φ β,j 2 = P β,k Φβ,k 2 L 2 (m β ) + Φβ,k , Φ β,0 2 + k-1 j=1 Φβ,j , Φ β,j Φβ,k , Φβ,j + Φβ,k , Φ β,j -Φβ,j , Φ β,j Φβ,j 2 ≤ P β,k Φβ,k 2 L 2 (m β ) + Φβ,k , Φ β,0 2 + 2 k-1 j=1 Φβ,k , Φβ,j 2 + 1 -Φβ,j , Φ β,j 2 = o β (1) + 2 k-1 j=1 1 -Φβ,j , Φ β,j 2 , (4.38)
where in the last line we used (4.32), (4.33) and (4.37).

In order to prove (ii), we now proceed by finite induction over k. The statement (ii) holds for k = 1 as the sum in the last line of (4.38) is actually empty. But then from (4.38) one immediately yields the inductive step.

ϕ dm β,k -ϕ Φ2 β,k dm β = Φ β,k , T ϕ Φ β,k -Φβ,k , T ϕ Φβ,k ≤ 4 ϕ C(M ) 1 -Φβ,k , Φ β,k 2 . ( 4
|Φ β,h Φ β,k |dm β ≤ E |Φ β,k | 2 dm β 1/2 + E c |Φ β,h | 2 dm β 1/2 , (4.43)
and both terms vanish as β → ∞.

I β (µ) ≥ β -2 dµ β 2 |∇V | 2 -β∆V -β 2 2 |∇V | 2 = 1 2 dµ |∇V | 2 -2 β ∆V ≥ 1 2 dµ R.
(5.1)

In particular, for any t > 0 and µ ∈ P(M ), it holds µ({R ≤ t}) ≥ 1 -2t -1 I β (µ). By A.2, R has compact sublevel sets, therefore any family (µ β ) ⊂ P(M ) such that (I β (µ β )) is uniformly bounded, is tight. And, by Prohorov's Theorem, precompact.

Proof of (2.13), Γ-liminf. Assume that µ β → µ. Then, replacing ψ by β ψ in (3.2), one obtains for each ψ ∈ Ω

I β (µ β ) ≥ dµ β ∇V • ψ -1 β ∇ • ψ -1 2 ψ 2 , (5.2)
which implies

lim β I β (µ β ) ≥ dµ ∇V • ψ -1 2 ψ 2 . ( 5.3)
Optimizing over ψ we get the result.

Proof of (2.13), Γ-limsup. We proceed in three steps.

Step1. First assume µ = δ x for some x ∈ M . Take U ≡ U x ∈ C 2 (M ; [0, +∞)) a smooth function such that (i) U (x) = distance(x, x) 2 in a neighborhood of x, and U is strictly positive elsewhere.

(ii) If M = R d , then U (x) ≥ |x| 2 + |∇V | 2 (x) + |∆V |(x) for x large enough and ∇ exp(-U ) ∈ L 2 (dx) (it's immediate to check that such a U exists). Let then µ β,x (dx) = 1 C β exp(-2β U (x))dx. ( 5.4)
By Laplace principle µ β,x → δ x. On the other hand an explicit calculation shows

I β (µ β,x ) = dµ β,x |∇V | 2 2 + 2|∇U | 2 - 1 β ∆V . (5.5) If M is compact, then lim β dµ β,x |∇U | 2 = |∇U | 2 (x) = 0, (5.6) lim β dµ β,x ∆V = ∆V (x), (5.7)
so that the second and third term in the l.h.s. of (5.5) vanish to get Step2.

lim β I β (µ β,x ) = lim β dµ β,x |∇V | 2 2 = |∇V | 2 (x) 2 = I(δ x). ( 5
If µ = N k=1 α k δ xk , then define µ β = N k=1 α k µ β,x k
, where µ β,x is defined as in (5.4). Clearly µ β → µ. On the other hand, since I β is convex and I is linear, again lim β I β (µ β ) ≤ I(µ).

Step3. By a standard diagonal argument, it's now enough to show that finite convex combinations of Dirac masses are I-dense in P(M ). This is trivial, since I is linear.

First order development by Γ-convergence

In this section we prove (2.14). The following remark is a restatement of (3.2). Remark 6.1. If µ = m β and h ∈ C 2 (M ) is such that h 2 = , then

I β (µ) = - 2 β 2 dm β (x)h(x)(L β h)(x). (6.1)
In particular if dµ = h 2 dm β and h is an L 2 (m β )-normalized eigenfunction of -1 β L β with eigenvalue λ, then

I β (µ) = 2 β λ. (6.2)
More generally, for P β,ε,λ defined as in (4.5) and for h ∈ Range(P β,ε,λ )

2 β (λ -ε) ≤ I β (µ) ≤ 2 β (λ + ε). (6.3)
Proof of (2.14): Γ-liminf inequality. Let µ β be a sequence converging to µ in P(M ), we need to show that lim

β βI β (µ β ) ≥ J(µ). (6.4)
From the Γ-liminf (2.13) proved in Section 5, the liminf in (6.4) equals +∞ unless I(µ) = 0, namely µ has the form

µ = z∈℘ α z δ z . (6.5)
Up to passing to a suitable subsequence (still labeled β hereafter), we can assume that the liminf is actually a limit, so that with no loss of generality I β (µ β ) < +∞ for all β large enough. Then by Proposition 3.3 we can assume that β = dµ β /dm β is in C 2 (M ) and bounded away from 0. In particular h

β := √ β ∈ C 2 (M ). Fix ε > 0, Λ large enough and such that Λ > max z∈℘ λ z,0 , (6.6)
and let β ≥ max (z,n) : λz,n≤Λ β ε,z,n . By Proposition 4.2 and using the notation therein introduced, recalling also (4.3)-(4.5), we can write

h β = λ≤Λ P β,ε,λ h β + P β,Λ h β , (6.7)
where

P β,Λ := 1 (Λ,+∞) (-1 β L β ) , P β,ε,λ h β = (z,n)∈S λ γ β,z,n Ψ β,z,n , γ β,z,n := h β , Ψ β,z,n ∈ [-1, 1].
By (6.1), Proposition 4.2 and the Markov spectral inequality

βI β (µ β ) = 2 λ≤Λ -1 β L β P β,ε,λ h β , P β,ε,λ h β m β + 2 -1 β L β P β,Λ h β , P β,Λ h β m β ≥ 2 λ≤Λ (λ -ε) P β,ε,λ h β 2 L 2 (m β ) + 2Λ P β,Λ h β 2 L 2 (m β ) = 2 λ≤Λ (z,n)∈S λ γ 2 β,z,n (λ z,n -ε) + 2Λ P β,Λ h β 2 L 2 (m β ) ≥ 2 z∈℘ α β,z,Λ (λ z,0 -ε) + 2Λ P β,Λ h β 2 L 2 (m β ) , (6.8)
where we used in the last line λ z,n ≥ λ z,0 and introduced

α β,z,Λ = n : λz,n≤Λ γ 2 β,z,n . (6.9)
As we take β → ∞, (6.8) holds for any ε > 0, in particular we gather

lim β P β,Λ h β 2 L 2 (m β ) ≤ lim β 1 2Λ βI β (µ β ) (6.10)
and, since

z∈℘ α z = 1 = h β 2 L 2 (m β ) = P β,Λ h β 2 L 2 (m β ) + z∈℘ α β,z,Λ , (6.11)
still by (6.8) and the fact that ε is arbitrary (once we take β → ∞)

lim β βI β (µ β ) ≥ lim β 2 z∈℘ α β,z,Λ λ z,0 + 2Λ z∈℘ α z -α β,z,Λ = 2 z∈℘ α z λ z,0 + 2 lim β z∈℘ (α z -α β,z,Λ )(Λ -λ z,0 ) = J(µ) + 2 lim β z∈℘ (α z -α β,z,Λ )(Λ -λ z,0 ), (6.12)
where in the last line we used (4.2). Since the sum above is finite, in view of (6.6), we conclude once we show that lim

β α β,z,Λ ≤ α z , ∀z ∈ ℘. (6.13)
To this aim, fix z ∈ ℘ and χ ∈ C b (M ), 0 ≤ χ ≤ 1, such that χ(z) = 1 and χ(z ) = 0 for z ∈ ℘, z = z. Then

α z = µ(χ) = lim β µ β (χ) = lim β m β (h 2 β χ) = lim β λ,λ ≤Λ (z,n)∈S λ (z ,n )∈S λ γ β,z,n γ β,z ,n χΨ β,z,n , Ψ β,z ,n L 2 (m β ) + lim β χP β,Λ h β , 2h β -P β,Λ h β L 2 (m β ) . (6.14)
Using (4.7) both for the on-and off-diagonal terms in the above sum and recalling (6.9) we obtain

α z ≥ lim β α β,z,Λ -2 χ C b P Λ,β h β L 2 (m β ) ≥ lim β α β,z,Λ -lim β 2βI β (µ β ) Λ , (6.15)
where in the last inequality we used (6.10). Since we can assume lim β βI β (µ β ) = C < +∞ with no loss of generality, and since α β,z,Λ is increasing in Λ, from (6. As we take the limsup β → +∞, (6.20) holds for any ε > 0, and thus lim β I β (µ β ) ≥ I(µ).

Exponential development by Γ-convergence

In this section we prove (2.15).

Proof of (2.15): Γ-liminf inequality. Fix some k = 1, . . . , n and let µ β be a sequence converging to µ. We need to show that lim β βe βW k I β (µ β ) ≥ J k (µ). (7.1) Up to passing to a suitable subsequence (still labeled β hereafter), we can assume that the liminf is actually a limit, so that with no loss of generality I β (µ β ) < +∞ for all β large enough. Then by Proposition 3.

Proposition 3 . 3 .

 33 Assume that Z := dx e -V (x) < +∞ and V ∈ C 1 (M ). Then I is the lower semicontinuous envelope of Ĩ : P(M ) → [0, +∞] defined as

. 17 )

 17 Moreover, from (4.14) and L 2 (dx) orthogonality of the Hermite functions (Θ z,n ) n it easily follows for ϕ ∈ C b (M ) lim β Θβ,z,n Θβ,z ,n ϕ dx = ϕ(z) if z = z and n = n 0 otherwise . (4.18) Finally define Ψβ,z,n = Θβ,z,n e βV 2

1 β

 1 .30) V (x k ) is the unique minimum of V (x) on B k , and χ(x k ) = 1. Therefore (i) holds. By Laplace principle lim β log Z β,k = -V (x k ), (4.31) lim β Φβ,j , Φβ,k = 0 j, k = 1, . . . , n, j = k, (4.32) lim β Φ β,0 , Φβ,k = 0 k = 1, . . . , n.

Proposition 4 . 7 .

 47 For k = 0, . . . , n let Φ β,k be a normalized eigenfunction corresponding to β,k and define the probability measure m β,k bydm β,k = Φ 2 β,k dm β . (4.39) Then m β,k converges to δ x k in P(M ) as β → ∞. In particular for h = k it holds lim β |Φ β,h Φ β,k |dm β = 0. (4.40) Proof. The statement (4.39) is immediate for k = 0, since we assumed x 0 to be the unique global minimizer of V . Let k ≥ 1, take Φβ,k as in Lemma 4.6 and fix a test function ϕ ∈ C b (M ). Now apply Remark 4.4 with H = L 2 (m β ), v = v = Φ β,k , u = u = Φβ,k , and T ≡ T ϕ the pointwise multiplication by ϕ. We get from (4.21)

5. Γ-convergence of I β Proposition 5 . 1 .

 51 Assume V ∈ C 2 (M ) and A.2. Then (I β ) is equicoercive. Proof. Let R(x) := inf β≥β0 |∇V | 2 (x) -2β ∆V (x). Taking ψ = β ∇V in the first line of (3.2) and recalling(3.6), one obtains for β ≥ β 0 and each µ ∈ P(M ),

  Proof of (2.14): Γ-limsup inequality. If J(µ) = +∞ there is nothing to prove. So we can assume µ = z∈℘ α z δ z . Recall the functions (Ψ β,z,n ) z∈℘,n∈N d introduced in Proposition 4.2 and define µ β,z ∈ P(M ) and the recovery sequence µ β ∈ P(M ) by dµ β,z := Ψ 2 β,z,0 dm β , (6.17)µ β := z∈℘ α z m β,z .(6.18)By (4.7), µ β,z → δ z in P(M ), so that µ β → µ.Using the notation of Proposition 4.2, for each ε > 0 and for each β ≥ max z∈℘ β ε,z,0 , by (6.3) and (4.2)I β (m β,z ) ≤ 2 β (λ z,0 + ε) = 1 β (ζ(z) + 2ε). (6.19) Now, since I β is convex by Remark 3.2 β I β (µ β ) ≤ β z∈℘ α z I(m β,z ) ≤ z∈℘ α z (ζ(z) + 2ε) = J(µ) + 2ε. (6.20)

2 h β 2 L 2 2 L 2 2 L 2

 2222222 [START_REF] Bovier | Metastability in reversible diffusion processes II. Precise asymptotics for small eigenvalues[END_REF] we can assume that β = dµ β /dm β is in C 2 (M ) and bounded away from 0. In particular h β := √ β ∈ C 2 (M ). As in the proof of Lemma 4.6 we denote by P β,k := 1 [ β,k+1 ,+∞) (-1 β L β ) the spectral projection of -1 β L β associated to the interval [ β,k+1 , +∞). Then, since h β ∈ L 2 (m β ), we can write h β = k j=0 γ β,j Φ β,j + P β,k h β , γ β,j = h β , Φ β,j ∈ [-1dm β (x) P β,k h β (x) (L β P β,k h β )(x) (m β ) ,(7.3)i.e.βe βW k I β (µ β ) ≥2γ 2 β,k e βW k β,k e βW k β,j + 2e βW k β,k+1 P β,k h β (m β ) .(7.4)From (4.28) and assumption A.5, the terms γ 2 β,j e βW k β,j with an index j ≤ k-1 vanish as β → ∞, since γ 2 β,j ≤ 1. On the other hand, still by (4.28) and A.5, the term e βW k β,k+1P β,k h β (m β ) diverges to +∞ unless P β,k h β L 2 (m β ) vanishes. Thus setting α := lim β γ 2β,k we gather again by (4.28)lim β βe βW k I β (µ β ) ≥ +∞ if lim β P β,k h β L 2 (m β ) > 0, αη k (x k ) otherwise. (7.5)In particular, if lim β P β,k h β L 2 (m β ) > 0, we proved the inequality (7.1) no matter what the limit µ of the µ β is. Now we claim that if lim β P β,k h β L 2 (m β ) = 0, then necessarily lim β µ β = µ = k j=0 α j δ xj for some α j , with α k = α, so that lim β βe βW k I β (µ β ) ≥ α k η k (x k ), namely (7.1). To prove this claim notice thatdµ β = (h β ) 2 dm β i,j=0 γ β,i γ β,j Φ β,i Φ β,j   dm β + (P β,k h β ) j Φ β,j + P β,k h β   dm β .(7.6)Since (P β,k h β ) L 2 (m β ) vanishesas β → ∞ and (4.40) holds, the last two term above also vanish (weakly) as measures as β → ∞. On the other hand, since µ β → µ, it follows from Proposition 4.7 that α j = lim β γ 2 β,j exists and µ = k j=0 α j δ xj . Proof of (2.15): Γ-limsup inequality. If J k (µ) = +∞ there is nothing to prove. So we can assume µ = k j=0 α j δ xj . For m β,j as in (4.39), define the recovery sequence µ β by µ β := k j=0 α j m β,j . (7.7) By Proposition 4.7, µ β → µ in P(M ). I β is convex, and by (6.2) I β (µ β ) ≤ 28) and hypotheses A.5 one finally gets lim β βe βW k I β (µ β ) ≤ 2 α k lim β e βW k β,k = α k η k (x k ) = J k (µ). (7.9)

  4.1. Quasimodes of the bounded spectrum. For z ∈ ℘ and n = (n 1 , . . . , n d ) ∈

	N d let
	λ z,n :=	ζ(z) 2	+	d k=1	n k |ξ k (z)|,	(4.2)

where ζ and the ξ k 's are defined in Subsection 2.4. The set (λ z,n ) n∈N d is nothing but the spectrum of a suitably rescaled quadratic approximation (i.e. of a harmonic oscillator) around the point z ∈ ℘ of the Schrödinger operator (1.8), see also the proof of Lemma 4.3 for this interpretation. Define also

  .24)In particular for β large enough the ( Ψβ,z,n ) (z,n)∈S λ are linearly independent, and span Range(P β,εΛ,λ ) by Theorem 4.1. Now construct (Ψ β,z,n ) (z,n)∈S λ by applying the Gram-Schmidt procedure on ( Ψβ,z,n ) (z,n)∈S λ . From (4.24) and Lemma 4.3-(ii) (with ϕ ≡ 1) it follows that lim

	0 otherwise	.	(4.25)
	Now apply Remark 4.4 with

β Ψβ,z,n , Ψ β,z ,n L 2 (m β ) = 1 if z = z and n = n

  The estimates of β,k+1 and Z β,k provided respectively in (4.28) and (4.31), and the choice (4.29) of δ thus imply

	lim β	P β,k Φβ,k	L 2 (m β )	= 0.	(4.37)

  → δ x k as β → ∞. Now (4.40) is an immediate consequence of this fact, since for E ⊂ M a neighborhood of x h not containing x k

	This vanishes by Lemma 4.6-(ii). Therefore by Lemma 4.6-(i)
	lim β	ϕ dm β,k = lim β	ϕ Φ2 β,k dm β (x) = ϕ(x k ),	(4.42)
	namely m β,k
				.41)

  .8) If M = R d , (5.6), (5.7) and (5.8) still hold if we restrict the integral to a (large) ball B containing x. Since e -2βU /C β vanishes pointwise out of B, it's enough to show that the integrands in (5.6), (5.7) and (5.8) are uniformly integrable. This is easily established by condition (ii) above.

This should not be confused with what is called large and moderate deviations in the literature, which is not related to metastability. Here the observables are fixed once and for all, and not rescaled. And yet they feature multiple large deviations principles.

Acknowledgements: GDG gratefully acknowledges the financial support of the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 614492, helpful discussions with Tony Lelièvre and Tom Hudson, and the warm hospitality at the University of Rome La Sapienza during his frequent visits.

This work has been carried out thanks to the support of the A*MIDEX project (n. ANR11IDEX000102) funded by the Investissements d'Avenir French Government program, managed by the French National Research Agency (ANR).

We are grateful to Lorenzo Bertini for introducing us to the problem.


Pascal Giorgi
email: pascal.giorgi@univ-perp.fr

Christophe Nègre
email: christophe.negre@univ-perp.fr

Équipe Dali

Lp2a

Thomas Plantard
email: thomaspl@uow.edu.au

Subquadratic Binary Field Multiplier in Double Polynomial System

Keywords: Binary Field Multiplication, Subquadratic Complexity, Double Polynomial System, Lagrange Representation, FFT, Montgomery reduction

We propose a new space efficient operator to multiply elements lying in a binary field F 2 k . Our approach is based on a novel system of representation called Double Polynomial System which set elements as a bivariate polynomials over F 2 . Thanks to this system of representation, we are able to use a Lagrange representation of the polynomials and then get a logarithmic time multiplier with a space complexity of O(k 1:31 ) improving previous best known method.

Introduction

Efficient hardware implementation of finite field arithmetic, and specifically of binary field F 2 k , is often required in cryptography and in coding theory [START_REF] Berlekamp | Bit-serial Reed-Solomon encoder[END_REF]. For example in elliptic curve cryptosystem [START_REF] Koblitz | Elliptic curve cryptosystems[END_REF][START_REF] Miller | Use of elliptic curves in cryptography[END_REF], the main operation is the scalar multiplication on the curve, which necessitates thousands of multiplications and additions over a finite field. Similarly, hundreds of multiplications over a binary field are required for the Diffie-Hellman Key exchange protocol [START_REF] Diffie | New directions in cryptography[END_REF].

Previously to this work, several architectures have already been proposed to efficiently implement the arithmetic in F 2 k . These architectures are mostly dedicated to the multiplication since this operation is extensively used and is often the most expensive. Each of them takes advantage of a special representation of the field. In particular, one of them uses polynomial basis or shifted polynomial basis [START_REF] Mastrovito | VLSI architectures for computations in Galois fields[END_REF][START_REF] Koc | Mastrovito Multiplier for All Trinomials[END_REF][START_REF] Fan | Fast bit-parallel GF(2 n ) multiplier for all trinomials[END_REF] while another uses normal basis [START_REF] Gao | Normal Basis Over Finite Field[END_REF][START_REF] Koc | Low Complexity Bit-Parallel Canonical and Normal Basis Multiplier For a Class of Finite Fields[END_REF]. The latter providing a really efficient squaring in the field since in this basis the squaring is just a cyclic shift of the coefficients.

In these representations the main approach to perform the multiplication consists to express the ope-ration as a matrix-vector product with binary entries. Parallel architectures are thus capable to perform this product within logarithmic time. However, these architectures still achieve a space complexity of k 2 . According to the recent improvements proposed in [START_REF] Fan | A new approach to subquadratic space complexity parallel multipliers for extended binary fields[END_REF], one can still perform the matrixvector product in logarithmic time but with a space complexity of k 1:56 or k 1:63 . This has been made possible thanks to structured matrices such as Toeplitz ones and a divide-and-conquer approach for the products.

In this paper we propose a new approach which reduces the exponent in the space complexity to 1:31 while keeping a logarithmic time complexity. First, we introduce a novel system of representation, the Double Polynomial System. In this representation, elements of F 2 k are polynomials in two variables

A(t;Y ) = ∑ n 1 i=0 a i (t)Y i
where a i (t) have degree strictly less than r.

Therefore, as in classical polynomial representation, the multiplication can be performed in two steps : a polynomial multiplication, and then a reduction phase to reduce the degrees in Y and in t.

The reduction in Y is simple due to the definition of DPS. The same is not true for the reduction in t.

Here, we use a Montgomery-like reduction approach in order to perform this reduction with few polynomial multiplications, this enabling us to easily use the Fast Fourier Transform. Therefore, our multiplier fully benefits from the FFT process which is highly parallelizable and provides a subquadratic space complexity.

Hence, we propose a binary field multiplier which has a delay of (16 log 3 (k) + 20)T X + 8T A and a space complexity of O(k 1:31 ), where T X and T A correspond respectively to the delay of one XOR gate and one AND gate.

Let us briefly give the outline of the paper. We first introduce the DPS representation for binary fields F 2 k (Section 2). We present the DPS multiplication in Section 3 and discuss the problem of finding a suitable polynomial to achieve our Montgomery-like coefficient reduction in Section 4. Then, we present in Section 5 a modified version of our multiplication introducing Lagrange basis. We recall in Section 6 some basic facts on the architecture design of a ternary FFT on which we rely for our multiplier. We finally conclude this paper by a detailed explanation of the complete architecture for our DPS-Lagrange multiplier and its complexity analysis and comparison (Section 7).

DPS Representation

A binary field F 2 k is generally constructed as the set of polynomials modulo an irreducible polynomial

P P F 2 [t] of degree k F 2 k = F 2 [t]=(P(t)) = fA(t) P F 2 [t] s.t. deg A(t) < kg
We introduce a novel binary field representation, the Double Polynomial System (DPS), inspired from AMNS number system of [START_REF] Bajard | Modular number systems : Beyong the mersenne family[END_REF].

Definition 1 (DPS representation

A(t) 0 t 2 t 3 + t 2 + t + 1 t 3 + t A B 0 (t + 1)Y 2 Y + 1 Y 2 + t + 1 A(t) 1 t 2 + 1 t 3 + t 2 + t t 3 + t + 1 A B 1 (t + 1)Y 2 + 1 Y Y 2 + t A(t) t t 2 + t t 3 + t 2 t 3 + 1 A B t Y 2 + Y + 1 Y + t Y 2 A(t) t + 1 t 2 + t + 1 t 3 + t 2 + 1 t 3 A B t + 1 Y 2 + Y Y + t + 1 Y 2 + 1
In particular, we can verify that if we evaluate

(t + 1)Y 2 + 1 in γ, we get (t + 1)γ 2 + 1 = (t + 1)(t 3 + t 2 + t) 2 + 1 = t 2 + 1 mod P, as expected. One can also see that deg Y ((t + 1)Y 2 + 1) = 2 < 3 = n and deg t ((t + 1)Y 2 + 1) = 1 < 2 = r.

¦

Remark 1. The DPS can be seen as a generalization of the polynomial representation of double extensions F 2 rn . Such extensions are usually construc- ted first as F 2 r = F 2 [t]=(P(t)) and then as F 2 rn = F 2 r [Y ]=(Y n λ)) with λ P F 2 r , see [START_REF] Guajardjo | Efficient algorithms for elliptic curve cryptosystems[END_REF]. However, this construction is not possible when the degree k of the field F 2 k is prime. DPS provides an alternative for double extension in this situation. Remark 2. As in classical polynomial representation, the addition in DPS is just a parallel bitwise XOR on the coefficients.

We proceed now by considering the problem of the multiplication of two elements expressed in a DPS. This can be done in two steps as described in Algorithm 1.

The first step of the algorithm consists of a classical polynomial multiplication modulo the binomial

E(Y ) = (Y n λ).
C = AB mod (Y n λ). 2. Coefficients reduction : R = RedCoe f f (C).
It is clear from the DPS system and from the multiplication modulo a binomial Y n λ that C has coefficients c i (t) with degree in t bounded by 2(r 2) + deg t λ. Therefore, these coefficients must be reduced to get the result of the multiplication expressed in the DPS representation.

Multiplication in DPS

A straightforward method for the reduction phase in t of Algorithm 1 is to perform an Euclidean division

C = Q ¢ M + R where deg t R < r. This reduction is only valid if M(t;Y ) is monic in t and satisfies M(t; γ) 0 mod P(t) with deg t (M) = r: (1)
Generally, one can easily compute a polynomial M satisfying equation ( 1), e.g. Section 4, but ensuring monicity is difficult.

In order to avoid monicity attached to a division strategy, we adapt the Montgomery trick (Montgomery, 1985) to our DPS system. The idea is to replace the Euclidean division by few multiplications and one exact division. This corresponds to annihilating the lower part of the c i (t) instead of the higher ones.

This method is given in Algorithm 2 assuming a polynomial M(t;Y ) satisfying M(t; γ) 0 mod P(t) is given.

Algorithm 2 DPS Multiplication.

Input : A; B P B = (P; γ; n; r; λ) with E = Y n λ Data : M such that M(γ) 0 mod P, a polynomial m P F 2 [t]
and

M H = M 1 mod (E; m) Output: R such that R(t; γ) = A(t; γ)B(t; γ)m 1 mod P begin C 2 A ¢B mod E; Q 2 C ¢M H mod (E; m); R 2 (C + Q ¢M mod E)=m;
end Example 2. We consider the field F 2 4 , with the DPS

B = (P = t 4 + t 3 + t 2 + t + 1; γ = t 3 + t 2 + t; n = 3; r = 2; λ = t).
In Table 2, we give an example of trace of DPS multiplication.

TAB. 2: DPS multiplication trace.

Operations Results

A tY 2 + tY B (t + 1)Y + t M tY 2 + Y + t + 1 M H (1 + t)Y 2 + (1 + t)Y + 1 m t 2 C tY 2 + t 2 Y + t 3 + t Q tY 2 Q ¢M (t 2 + t)Y 2 + t 3 Y + t 2 C + Q ¢M t 2 Y 2 + (t 3 + t 2 )Y + t 3 R Y 2 + (t + 1)Y + t
We can check that R(t; γ) t 2 + t mod P is equal to A(t; γ)B(t; γ)t 2 mod P.

¦

Lemma 1. Algorithm 2 is correct.

Proof. We need to demonstrate that the output R of the algorithm satisfies the following equation R(t; γ) = A(t; γ)B(t; γ)m 1 mod P:

(2)

From the definition 1 of DPS representation, we know that E(γ) 0 mod P. Thus, we have C(t; γ) A(t; γ)B(t; γ) mod P: By definition of M, we have M(t; γ) 0 mod P and consequently

C(t; γ) + Q(t; γ)M(t; γ) C(t; γ) A(t; γ)B(t; γ) mod P
We now need to prove that the division by m is exact. This is equivalent to prove the following equi-

valence C + Q ¢M mod E 0 mod m. By definition, we have Q = C ¢ M H mod E and M H = M 1 mod (E; m). We consider R H = C + Q¢M mod (E; m), then the following equivalences hold R H C + C ¢( M 1 ¢M) mod (E; m) (C C) mod (E; m) 0 mod (E; m):
Thus, division by m is exact. Hence, the algorithm is correct since an exact division (the division by m) is equal to the multiplication by an inverse modulo P.

At this level, we know that the resulting polynomial R of the previous algorithm satisfies the equation R(t; γ) = A(t; γ)B(t; γ)m 1 mod P but we do not know whether it is expressed in the DPS, i.e., if the coefficients of R have degree in t smaller than r. This is the goal of the following theorem.

Construction of the polynomial M

The result of this section uses mathematical structures involving module over the polynomial ring F 2 [t] in order to prove existence of a suitable polynomial M. The remaining of the paper is independent from this section and readers who are not familiar with such mathematical structure can skip this section without misunderstanding.

Our goal is to construct a polynomial M such that M(t; γ) 0 mod P and deg t M is small. This polynomial belongs to the set

M = fA(t;Y)) P F 2 [t;Y ] with deg Y A < ng:
The set M has a natural structure of F 2 [t] module. Recall that an F 2 [t]-module M is an (additive) abelian group, with a scalar multiplication over F 2 [t] :

F 2 [t] ¢M 3 M :
In order to calculate the element M with low degree in t, we will use a sub-module M H of M spanned by the following linearly independent vectors. Ω = 0 B B B B B @ P 0 0 : : : 0 γ 1 0 : : : 0 γ 2 0 1 : : : 0 . . . . . . . . .

γ n 1 0 0 : : : 1 1 C C C C C A 2 P 2 Y γ 2 Y 2 γ 2 . . . 2 Y n 1 γ n 1
Each of the polynomials V (t;Y ) defined by the rows of Ω satisfy V (t; γ) 0; and any F 2 [t]-linear combination of these polynomials satisfies also this property. Therefore, one way to construct M consists to compute a minimal basis of M H and define M as the basis element with the smaller degree in t. The notion of minimality is related to the degree in t of the basis elements.

According to polynomial matrix properties, one can find a minimal basis of Ω by computing its matrix reduced form called the Popov form [START_REF] Mulders | On lattice reduction for polynomial matrices[END_REF]. In particular, the properties of the Popov form (Villard, 1996, §1.2) tell us that it exist a minimal basis ( f 1 ; f 2 ; : : : ; f n ) of M H which satisfies the following degree properties : Beside the fact that the calculation of M is only needed once at the construction of the DPS representation, one would need to efficiently compute such polynomial. This can be achieve within a complexity of O(n 3 k 2 ) binary operations with Algorithm WeakPo-povForm of [START_REF] Mulders | On lattice reduction for polynomial matrices[END_REF] or with an asymptotic complexity of O(n 3 k log k) binary operations with Algorithm ColumnReduction of [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF].

n ∑ i=1 deg t f i = deg t (det(Ω)) (4)

Lagrange representation

Let R a ring, and R [Y ] the polynomial ring over R . The Lagrange representation of a polynomial of degree n 1 in R [Y ] is given by its values at n distinct points. For us, these n points will be the roots of a polynomial

E = ∏ n i=1 (Y α i ) P R [Y ].
From an arithmetic point of view, this is related to the Chinese Remainder Theorem which asserts that the following application is an isomorphism

R [Y ]=(E(Y )) g 3 n ∏ i=1 R [Y ]=(Y α i ) (7) A U 3 (A mod (Y α i )) iPf1;:::;ng :
The computation of A mod (Y α i ) is simply the computation of A(α i ). In other words, the image of A(Y ) by the isomorphism ( 7) is nothing else than the multi-points evaluation of A at the roots of E. This fact motivates the following Lagrange representation of the polynomials.

Definition 2 (Lagrange representation). Let

A P R [Y ] with deg A < n,

E(Y

) = r ∏ i=1 (Y α i ) mod m If a i = A(α i ) for 1 i n, the Lagrange repre- sentation (LR) of A(Y ) is defined by LR(A(Y )) = (a 1 ; : : : ; a n ):
Lagrange representation is advantageous to perform operations modulo E : this is a consequence of the Chinese Remainder Theorem. Specifically the arithmetic modulo E in classical polynomial representation can be costly if E has a high degree. In LR representation this arithmetic is decomposed into n independent arithmetic units, each does arithmetic modulo a very simple polynomial (X α i ). Furthermore, arithmetic modulo (X α i ) is the arithmetic in R since the product of two zero degree polynomials is just the product of the two constant coefficients.

Multiplication algorithm

Let us go back to the Algorithm 2 and see how to use Lagrange representation to perform polynomial arithmetic in each step. The first two steps can be done in Lagrange representation modulo m 1 (t) such that E split modulo m 1 (t) :

E = n ∏ i=1 (Y α i ) mod m 1 (t);
The third step must be done modulo a second polynomial m 2 (t), which also splits E,

E = n ∏ i=1 (Y α H i ) mod m 2 (t);
since the division by m 1 cannot be performed modulo the polynomial m 1 (t).

We then need to represent the polynomials A and B in Algorithm 2 with both their Lagrange representations modulo m 1 (t) and m 2 (t).

Notation 1. We will use in the sequel the following notation. For a polynomial A of degree n 1 in Y we will denote -A the Lagrange representation in α i modulo m 1 (t)

-A the Lagrange representation in α H i modulo m 2 (t).
Hence, we can do the following modifications to the Algorithm 2 :

Algorithm 3 DPS-LR Multiplication. Input : A; A; B; B Data : M such that M(t; γ) 0 mod P, M H such that M H = M 1 (mod E; m 1 ). Output: R; R such that R P B and R(t; γ) = A(t; γ)B(t; γ)m 1 1 mod P(t) begin Q 2 A ¢B¢M H ; Q 2 Convert m 1 3m 2 (Q)); R 2 (A ¢B)+Q¢M)¢m 1 1 ; R 2 Convert m 2 3m 1 (R); end
The operations to compute Q and R are performed in Lagrange representation and then can be easily parallelized. It consists of n independent multiplications

in F 2 [t]=(m 1 (t)) and F 2 [t]=(m 2 (t)).
The major drawback of this algorithm is the conversions between Lagrange representations modulo m 1 and m 2 . It is necessary to perform these operations efficiently in order to get a multiplier yielding our announced space complexity.

Lagrange representations conversion

In order to provide an efficient implementation of conversions between Lagrange representations modulo m 1 and m 2 , we rely on the binomial form of

E = Y n λ. Indeed, if µ 1 = α 1 is a root of E modulo m 1
then all others roots can be written α j = µ 1 ω i 1 mod m 1 where ω 1 is a n-th primitive root of unity in F 2 [t]=(m 1 ). This property comes from the fact that (α j =µ 1 ) n = 1 mod m 1 and thus there exists an integer i such that α j =µ 1 = ω i 1 mod m 1 . This is still true modulo m 2 .

Thus, the multi-point evaluation of the polynomial A(Y ) in α i modulo m 1 can be done as follow :

1. set e A(Y ) = A(µ 1 1 Y ) = n 1 ∑ i=0 a i µ i 1 Y i 2. compute A = DFT m 1 ( e A; n; ω 1 ),
where DFT m 1 ( e A; n; ω 1 ) is the evaluations of the polynomial e

A in the n-th roots of unity ω i 1 . Similarly the Lagrange interpolation which compute A(Y ) from A can be done by reversing the previous process.

By gluing together this two processes we get the following algorithm to perform conversion between Lagrange representations.

Algorithm 4 Convert m 1 3m 2 . Input : A Output: A e A(Y ) 2 DFT 1 m 1 (A; n; ω 1 ) ; A(Y ) 2 e A(µ 1 1 Y ) mod m 1 ; e A(Y ) 2 A(µ 2 Y ) mod m 2 ; A 2 DFT m 2 ( e A(Y ); n; ω 2 );
As a consequence, the conversion has a the cost of two Discrete Fourier Transforms. This can be done efficiently by using FFT algorithm (Gathen and Gerhard, 1999, §8.2). Moreover, such FFT computations can be done even faster when special moduli m i (t) such that m i (t) = t 2¢3 s + t 3 s + 1 are used [START_REF] Schönhage | Schnelle multiplikation von polynomen über körpern der charakteristik 2[END_REF]. These moduli can be seen as a polynomial analog of Fermat numbers.

Architecture for FFT computation

We present an architecture to perform the FFT calculation of a polynomial A(Y ) P R [Y ] of degree n 1, keeping in mind our targeted Lagrange conversion algorithm. We consider the ring R

= F 2 [t]=(m(t))
where m(t) = t 2n=3 + t n=3 + 1 and n = 3 s . Note that the FFT process needs to be performed using ternary method since binary one is not feasible over characteristic 2 rings [START_REF] Schönhage | Schnelle multiplikation von polynomen über körpern der charakteristik 2[END_REF].

Let us denote ω a primitive n-th root of unity modulo m(t) and θ = ω n=3 a 3rd root of unity. The ternary FFT process is based on the following three-way splitting of A

A 1 = ∑ n=3 1 j=0 a 3 j Y 3 j ; A 2 = ∑ n=3 1 j=0 a 3 j+1 Y 3 j ; A 3 = ∑ n=3 1 i=0 a 3 j+2 Y 3 j ; such that A = A 1 + YA 2 + Y 2 A 3 . Let b A[i] = A(ω i ) be the i-th coefficient of DFT m (A; n; ω).
Let us also denote by Â1 [i]; Â2 [i] and Â3 [i] the coefficients of the DFT of order n=3 of respectively A 1 ; A 2 and A 3 .

The following relations can be obtained by evaluating

A = A 1 + YA 2 + Y 2 A 3 in ω i ; ω i+n=3 and ω i+2n=3 : Â[i] = Â1 [i] + ω i Â2 [i] + ω 2i Â3 [i]; Â[i + n=3] = Â1 [i] + θω i Â2 [i] + θ 2 ω 2i Â3 [i]; (8) Â[i + 2n=3] = Â1 [i] + θ 2 ω i Â2 [i] + θω 2i Â3 [i]:
This operation is frequently called the butterfly operation. It can be performed efficiently if we compute modulo m(t)(t n=3 + 1) = t n + 1 instead of m(t).

Indeed, in this case ω = t and a multiplication a(t) ¢ ω i modulo t n + 1 is a simple cyclic shift. The butterfly circuit (Figure 1) is a consequence of this remark and the relations given in (8). In Figure 1, the ) blocks refer to a simple shift operations by the given value and the L blocks refer to XOR operator. When no value is given, then shift operation is not performed.

2(i + n=3) i + n=3 2i + n=3 i + 2n=3 i 2i Â[i + n 3 ] Â[i + 2n 3 ] Â[i] Â1 [i]
Â2 [i] Â3 [i] FIG. 1: Ternary butterfly operator.

Within the FFT, the computations of Â1 ; Â2 and Â3 are done in the same way. These polynomials are split in three parts and butterfly operations are applied again. This process is done recursively until constant polynomial are reached. If we entirely develop this recursive process we obtain the schematized architecture in Figure 2. Let us now evaluate the complexity of this architecture. It is composed of log 3 (n) stages where each stage consists of n=3 butterfly operations. Each of these butterfly operations requires 6n XOR gates, and has a delay of 2T X , where T X is the delay of one XOR gate.

Consequently, this architecture has a space complexity of

S(FFT

m(t) ) = (2n log 3 (n) + n) XOR (9)
and a delay of

D(FFT m(t) ) = (2 log 3 (n) + 1)T X : (10)

Architecture and Complexity

We now present an hardware architecture associated to Algorithm 3 in the special case where m 1 and m 2 are chosen as

m 1 = t 2n + t n + 1 and m 2 = t 2n=3 + t n=3 + 1:
This choice enables us to use the FFT circuit presented in the previous section. The architecture of our binary field multiplier is given in Figure 3. It is constituted of FFT blocks and multipliers modulo m 1 (t) and m 2 (t).

TAB. 3: Complexity of multipliers modulo m 1 and m 2 .

Mul m1

Space Time #AND 3n log 3 (6) 1 #XOR 72 5 n log 3 (6) 9n 7=5 3 log 3

(n) + 3 Mul m2 Space Time #AND 1 2 n log 3 (6) 1 #XOR 36 15 n log 3 (6) n=5 + n 1 3 log 3 (n)
These multipliers are referenced by blocks Mul m 1 and Mul m 2 in our architecture. Because of the special form of m 1 (t) and m 2 (t) we can use the multiplier of [START_REF] Fan | A new approach to subquadratic space complexity parallel multipliers for extended binary fields[END_REF] to perform this operation. Therefore, the complexity (cf. Table 3) of these blocks are easily deduced from [START_REF] Fan | A new approach to subquadratic space complexity parallel multipliers for extended binary fields[END_REF], Table 1).

The FFT blocks are designed using ternary method presented in previous section. Therefore, their complexity are those given in ( 9) and (10).

1 µ 1 1 µ (n 1) 1 µ n 1 1 µ (n 1) 2 a 1 m H 1 b 1 b n 1 a n 1 FFT m1 b 1 a 1 m 1 b 1 a 0 m 0 b n 1 m n 1 a n 1 µ 2 µ n 1 2 1 1 1 µ 1 µ 1 2 Mul m1 Mul m1 Mul m1 Mul m1 m H 0 b 0 a 0 Mulm 2 Mulm 2 Mulm 2 Mulm 2 Mulm 2 Mulm 2 Mulm 2 Mulm 2 Mulm 2 Mulm 2 m 1 1 Mul m1 Mul m1 Mul m1 Mul m1 Mul m1 Mul m1 Mul m2 A ¢B¢M H convert m23m1 FFT m1 m H n 1 r 0 r 1 r n 1 FFT m2 FFT m2 Mul m2 Mul m2 Mul m2 Mul m1 Mul m2 Mul m2 Mul m1 r 0 r 1 r n 1 m 1 1 m 1 1 (A ¢B+Q¢M)m 1 1 convert m13m2 FIG. 3: DPS-Lagrange Multiplier.
The complexity of our multiplier can be evaluated with respect to the numbers of each blocks and their corresponding space complexity denoted S, and time complexity denoted D. For the space complexity this gives 4nS (Mul m 1 ) + 5nS (Mul m 2 ) + 2S (F FT m 1 ) + 2S (F FT m 2 ) + 2n 2 =3 XOR: Similarly, the critical path of this architecture gives the delay 4D(Mul m 1 ) + 4D(Mul m 2 ) + 2D(FFT m 1 ) + 2D(FFT m 2 ) + T X : Using these expressions, ( 9),(10) and FH £ = [START_REF] Fan | A new approach to subquadratic space complexity parallel multipliers for extended binary fields[END_REF] ; can compute the complexity with respect to the number of XOR and AND gates and their corresponding delay T X and T A .

Let r be the degree in t of the coefficients in the DPS representation then deg t (m 2 ) must satisfy deg t (m 2 ) ! r. Therefore, this implies that k r ¢n = 2n 2 =3 and thus leads to use n % p k, where k is the degree of the field F 2 k .

Finally, we obtain the complexity of the DPS-Lagrange multiplier stated in Table 4. We also give in this table the complexity of the best known method, regarding space and time complexity, to perform binary field multiplication.

One can remark that our approach decrease the space complexity from k 1:58 to k 1:31 , while it is slower by a factor roughly equals to 5:3.

For the sake of simplicity, we have presented our algorithm with special moduli m 1 and m 2 . As a consequence, we have an increase in the degree of m 1 which has an influence on the final complexity. Another approach would be to use only m 2 as a special Fermat polynomial and takes m 1 such that F 2 [t]=(m 1 ) and F 2 [t]=(m 2 ) are isomorphic. Therefore, operations modulo m 1 could be performed modulo m 2 through isomorphism transformations, and then decrease the space complexity by a factor at least 3.

Conclusion

In this paper we have presented a novel algorithm to perform multiplication in binary field, using a Double Polynomial System of representation. This system enables the use of Fast Fourier Transform in the multiplication according to Lagrange representation. The resulting multiplier still achieve a logarithmic time complexity, but asymptotically improve the space complexity from O(k 1:58 ) to O(k 1:31 ),

Our method is a first approach to reduce the space complexity of binary field multiplier. In particular, some optimizations can be done to reduce the constant factors in the complexity. For example, lot of multiplications by a constant are counted as full multiplication in the current complexity evaluation.

Furthermore, one can also reduce the exponent in the space complexity by replacing Fan and Hasan multipliers with a quasi-linear approach (e.g. Schönhage's technique [START_REF] Schönhage | Schnelle multiplikation von polynomen über körpern der charakteristik 2[END_REF]).

Theorem 1 .

 1 Let B = (P; γ; n; r; λ) a Double Polynomial System, M be a polynomial of B such that M(γ) 0 mod P and σ = deg t (M). Let A; B be two elements expressed in the DPS B. If r and the polynomial m satisfy r > σ + deg t (λ) and deg t (m) > deg t (λ) + r (3) then the polynomial R output by the Algorithm 2 is expressed in the DPS B.

  and α 1 ; : : : ; α n be the n distinct roots of a polynomial E(Y ).

  Proof. From the Definition 1, the polynomial R belongs to the DPS B = (P; γ; n; r; λ) if deg Y R < n and if deg t (R) < r. The fact that deg Y R < n is easy to see since all the computation in the Algorithm 2 are done modulo E = Y n λ. Hence, we have only to prove that deg t R < r. Since by definition deg t A; deg t B < r we have the following inequalities deg t R = deg t ((A ¢B+Q¢M) mod E)=m) max(deg t A + deg t B; deg t Q + deg t M) + deg t λ deg t m max(2r; σ + deg t m) + deg t λ deg t m: According to our hypothesis in the equation (3), we have both 2r + deg t λ deg t m < r and σ + deg t m + deg t λ deg t m < r. Hence, we get deg t (R) < r as required.

Table 3 , we

 3

  TAB. 4: Complexity comparison.

			Space Complexity	Time Complexity
	Method	# AND	# XOR	T A	T X
	This paper FH £ binary FH £ ternary	14:5k 1:31 69:6k 1:31 31k + k 0:5 (8 log 3 (k) + 39) 8 16 log 3 (k) + 20 k 1:58 5:5k 1:58 5k 0:5 1 2 log 2 (k) + 1 k 1:63 4:8k 1:63 4k 0:8 1 3 log 3 (k) + 1

DPS multiplication in Lagrange representationIn this section, we present a version of Algorithm 2 using a Lagrange representation of the DPS elements.


Keywords: Wavelet transforms, electrocardiograph, microwave system, cardiopulmonary signals

This paper presents different signal processing techniques used in separating cardiopulmonary signals measured using a microwave Doppler radar. For different sides from the person under test, the microwave system is tested for several operational frequencies simultaneously with a PC-based electrocardiograph. Wavelet transforms are used in order to separate heartbeat signal from the cardiopulmonary signals. High accuracy is obtained in terms of heartbeat rate for the whole measurements.

I. INTRODUCTION

he demand for wireless monitoring has increased recently for many applications including life sign detection for people under rubble or snow, newly born infant or burn victims, where electrocardiograph (ECG) with affixed electrodes can be perturbing. Thus, the utility of the microwave Doppler radar has increased as used for touch-less monitoring. According to the Doppler theory, a subject with a quasiperiodic movement will reflect the transmitted continuous wave signal with its modulated by the subject displacement [START_REF] Lin | Microwave sensing of physiological movement and volume change: A review[END_REF]. When transmitted to a person's chest, the signal will be reflected with the phase modulated according to the chest displacement. This latter is due to both, respiration and heart beating. As the heartbeat and the respiration signals are laid together, signal processing techniques are required to separate the heartbeat signal from the cardiopulmonary signal in order to extract the heartbeat rate (HR).

T

As the chest displacement due to respiration (4 -12 mm) is much higher than that due to heartbeat (0.2 -0.5 mm), the respiration rate can be determined without filtering while determining the heartbeat rate needs a processing technique [START_REF] Ramachandran | Three-dimensional reconstruction of cardiac displacement patterns on the chest wall during the p, qrs, and tsegments of the ECG by laser speckle interferometry[END_REF]- [START_REF] Kondo | Laser monitoring of chest wall displacement[END_REF]. This can be established due to the fact that the heartbeat rate varies between 50 and 90 beats per minute (0.83 to 1.5 Hz) while the respiration rate varies between 9 and 24 breaths per minute (0.15 to 0.4 Hz) [START_REF] Gorgas | Vital Signs and Patient Monitoring Techniques[END_REF]. This allows the extraction of an average heartbeat rate over a specific window of time. On the other hand, tracking the variations of the heartbeat rate over time requires signals separation and a peak-finding technique. Previous works tend to detect the heartbeat rate either by measuring separately the heartbeat and respiration signal, or by using classical filters, then by applying Fourier transform to the filtered signal and considering the peak to correspond to the HR [START_REF] Chen | Microwave lifedetection systems for searching human subjects under earthquake rubble or behind barrier[END_REF]- [START_REF] Lohman | A digital signal processor for Doppler radar sensing of vital signs[END_REF]. Other studies show the possibility of estimating the HR based on the periodicity of the autocorrelation function [START_REF] Lubecke | Wireless LAN PC card sensing of vital signs[END_REF] or upon applying FFT to the autocorrelated signal to obtain the heartbeat rate [START_REF] Li | Experiment and spectral analysis of a lowpower ka-band heartbeat detector measuring from four sides of a human body[END_REF], [START_REF] Li | Robust overnight monitoring of human vital signs by a noncontact respiration and heartbeat detector[END_REF].  Recently, another processing techniques are used for cardiorespiratory separation: field-programmable gate arrays (FPGAs) is used to process either time-or frequency-domain signals in human sensing radar applications applied for CW and UWB radars [START_REF] Wang | CW and Pulse-Doppler Radar Processing Based on FPGA for Human Sensing Applications[END_REF] as well as the use of compact quadrature Doppler radar sensor where Continuous-wavelet filter and ensemble empirical mode decomposition (EEMD) based algorithms are applied for cardiorespiratory signal to separate the cardiac and respiratory signals [START_REF] Hu | Noncontact Accurate Measurement of Cardiopulmonary Activity Using a Compact Quadrature Doppler Radar Sensor[END_REF]. These works show the possibility of detecting the HR for person when holding the breath as well as when breathing normally. However, they lack providing information about the heart rate variation in time as only an average value of the HR is extracted. The aim of this work is to show the possibility of detecting the cardiopulmonary signals of a person from four different sides: front, back, left, and right. The proposed system is tested at different operational frequencies: 2.4, 5.8, and 10 GHz. Simultaneously to the microwave system, PC-based ECG is used as a reference signal in order to validated the HR extraction. Obtained signals are processed using several wavelet transforms in order to extract the heartbeat signal along the time-axis. The rest of this paper is organized as follows: section II presents the proposed system and the measurement setup. Section III discusses the signal processing techniques used. Section IV shows and discusses the obtained results. Section V concludes the work.

II.PROPOSED SYSTEM AND MEASUREMENT SETUP

A. Proposed System

The contact-less detection of the cardiopulmonary activity is based on the reflected signal off the person's chest. The RF signal penetrates clothes with minimal reflection, and has a much higher reflection at the air/skin interface. This depends on the frequency of the transmitted signal. The signal reflection at the air/skin interface decreases as the frequency decreases, and has more significant reflections from clothing or bedding, and vice-versa. The amount of the phase variation is indirectly proportional to the wavelength of the carrier frequency. Thus, the signal-to-noise ratio is directly proportional as well to the operational frequency (f = C/λ). Higher frequency results in a shorter the wavelength, and then in a greater phase variation. Another feature related to the operational frequency is the size of the antenna. As the frequency increases, the same antenna gain can be achieved with a smaller antenna.

The proposed microwave system is based on using a Vector Network Analyzer (VNA), and two horn antennas. In addition to the installation simplicity, the VNA provide the ability of tuning both frequency and power of the transmitted signal. Moreover, it measures many parameters and especially the Advanced Signal Processing Techniques for Microwave Cardiopulmonary Signals Separation D. Obeid, S. Samad, S. Sadek, G. Zaharia, and G. El Zein time variation of the phase and module of the S parameters of the system under test. The phase of the S21 parameter corresponds to the difference between the phase of the received signal and the phase of the transmitted signal. Simultaneously with microwave system, a PC-based ECG is used in order to extract the electrocardiogram signal. This signal is used as a reference signal to which the processed microwave signal is compared. More details about the system setup can be found in [START_REF] Obeid | Position-Free Vital Sign Monitoring: Measurements and Processing[END_REF], [START_REF] Obeid | Microwave Doppler radar for heartbeat detection vs electrocardiogram[END_REF].

B. Measurement Setup

Measurements are performed on healthy subject of 54 years old, sitting at a distance of 1 meter from the antennas. As the life sign detection for patients and people under rubble requires heart activity detection regardless their positions with respect to the system, measurements are performed at the four different sides from the person under test: front, back, left, and right. As the reflection at the air/skin interface changes when the operational frequency changes, a comparative study is provided in terms of operational frequency where several frequencies are tested using the proposed system. The operational frequencies should cover as much radar band designations as possible, as well as ISM bands. Thus, the operational frequencies tested in this work are: 2.4 GHz (ISM S-band), 5.8 GHz (ISM C-band), and 10 GHz (X-band). From each side of the person under test, the system is tested at several operational frequencies: 2.4, 5.8, and 10 GHz. The total output power of the transmitted signal is 0 dBm. Each measurement lasts 30 seconds during which the subject breathes normally. Both the microwave and the ECG measurements start simultaneously. The phase variation due to the chest displacement is stored in the VNA and sent to the personal computer for processing. The processed signals are compared to the ECG signal. The ECG provides the electrical activity of the heart, however, only the RR interval is used as reference information for the comparison with the VNA filtered signal.

III.SIGNAL PROCESSING TECHNIQUES

As the phase variations of S21 due to breathing are higher than those caused by the heart beating, processing techniques are necessary to extract the heartbeat signal from the obtained cardiopulmonary signal. Previous work tends to apply the Fourier transform to the obtained cardiopulmonary signal in order to extract the heartbeat rate. This provides a mean value of the HR over a specific window of time; thus, it lacks providing the more important information about the variation of the HR in time, and cannot be established in real-time as it requires a window of long-duration. In order sto overcome these problems, the Discrete Wavelet Transform (DWT) is applied for the extraction the heartbeat signal.

The DWT ( W j ,k ) of a signal f(t) is given by the scalar product of f(t) with the scaling function (i.e. the wavelet basis function ϕ (t) which is scaled and shifted:

( W j ,k ) (1)
where the basis function is given by:

ϕ j ,k (t )=2 -j 2 ϕ (2 -j t-k ) (2)
where j is the j th decomposition level or step and k is the k th wavelet coefficient at the j th level [START_REF] Blackburn | Retrieval of chlorophyll concentration from leaf reflectance spectra using wavelet analysis[END_REF]. DWT is computed by successive low-and high-pass filtering of the discrete timedomain signal. The main advantage of the DWT determination is that the examination of the signal at different frequency bands with different resolutions by decomposing the signal into Approximation coefficients (A) and Detailed information (D). Hence, this algorithm gives precise analysis of frequency domain at low frequency and time domain at high frequency.

The "detail" Dn contains frequencies between fs/2 n and fs/2 n+1 . At rest, the HR varies between 60 and 120 beats per minute and the corresponding frequency is located between 1 and 2 Hz. The actual sampling frequency used in the VNA is 666.7 Hz, thus, there is no decomposition that provides the signal with frequency components between 1 and 2 Hz. A resampling is required in order to convert the sampling frequency from 666.7 Hz to 512 Hz which allows the extraction of the desired signal. This lets the 1-2 Hz components to be included in the 8 th decomposition of the wavelet. Once the wavelet decomposition is extracted, the signal is reconstructed in time-domain. A peak-detection method is applied to the time-domain reconstructed signal in order to detect the peaks that correspond to the beats and to extract the heartbeat rate.

Several wavelet families are tested in Matlab for the measurement performed at different sides from the person under test and for all the operational frequencies. The wavelet families include Bio2.4, Rbio2.3, Sym5, Db5, Coif3, and Dmey.

As not all the wavelet families fulfill the properties of the desired signal, the wavelet family has to be chosen as close as possible to the analyzed signal to give a better reconstruction with fewer decomposition levels [START_REF] Mohd Tumari | Selection of a Suitable Wavelet for Cognitive Memory Using Electroencephalograph Signal[END_REF], [START_REF] Megahed | Selection of a Suitable Mother Wavelet for Analyzing Power System Fault Transients[END_REF]. The choice of a suitable wavelet is based on the perfect reconstruction of the power system signal. The error between the original signal and the reconstructed signal should be the smallest for perfect reconstruction. The Mean Square Error (MSE) is calculated and used as indicator to measure the error between the original signal and the reconstructed signal. The MSE is computer with the following relation [START_REF] Mohd Tumari | Selection of a Suitable Wavelet for Cognitive Memory Using Electroencephalograph Signal[END_REF], [START_REF] Megahed | Selection of a Suitable Mother Wavelet for Analyzing Power System Fault Transients[END_REF]:

MSE=‖S-S reconstructed ‖= √ ∑ n= 0 N-1 |x (n )-x ( n) | 2 (3)
In order to compare the ECG signal to the reconstructed filtered cardiopulmonary signal, the peak detection method is applied to both signals. This allows comparing the beat-tobeat interval. The mean error is calculated for each measurement and for the different wavelet decompositions. The mean of the MSE of all available signals is calculated for each wavelet type according to the following equation: p×q where p indicates the positions and q corresponds to the emitted frequency. Four positions are considered: p=1 for the front side position, p=2 for the back side position, p=3 for the right side position, and p=4 for the left side position. Three emitted frequencies are considered during measurements: q=1 for 2.4 GHz, q=2 for 5.8 GHz, and q=3 for 10 GHz.

IV. RESULTS AND DISCUSSION

This section presents the results obtained upon measuring the cardiopulmonary signals using a microwave touch-less technique and for different sides from the person under test. Fig. 1 shows the cardiopulmonary signal detected from the four sides of the subject when using the 5.8 GHz operational frequency: a) the S21 variation due to breathing is shown when measured from the front side a), from the back side b), from the left side c), and from the right side d).

The measurement duration is 30 seconds during which the subject is requested to breathe normally. As seen in the figure, the time-domain respiration signal is clear when measured from the front side.

The signal separation is achieved upon applying several wavelet transforms families. Fig. 2 shows the D8 level extracted from the cardiopulmonary signal measured at 5.8 GHz from the front side upon using Coiflet wavelet (Coif3) compared to the corresponding ECG signal. where N is the number of peaks and dk is the duration between two consecutive peaks. The peaks that correspond to beats are detected using a peak-detection algorithm. The absolute value of the relative error of the HR is calculated as: Error= 100 * ( HR ECG -HR VNA ) HR ECG [START_REF] Chen | Microwave lifedetection systems for searching human subjects under earthquake rubble or behind barrier[END_REF] The following table shows the relative error calculated for the different wavelet transforms to each of the performed measurement. As shown in the table, most of the wavelet transforms show the ability to extract the heartbeat signal from the cardiopulmonary signal when using different operation frequencies and from all sides of the person under test. The lowest error obtained is when using the wavelet type Bior2.4.

Concerning the positions, it can be noticed that the right position results in the highest relative error (>10%). This is due to the fact the far side from the heart position is less affected by heart beating. All tested operational frequencies, 2.4, 5.8, and 10 GHz, show the possibility to track the cardiopulmonary signals from all sides.

V. CONCLUSION This paper presents the results of detecting cardiopulmonary signals from different sides of a person while breathing normally. Several wavelet transform families are used for signal separation and heartbeat rate extraction. Bior2.4 shows the highest accurate HR when compared to the ECG reference signal. Measurements are performed for several operational frequencies, 2.4, 5.8, and 10 GHz, for a distance of 1 meter from the subject, simultaneously with an electrocardiograph.

Future work will focus on tracking the heart activity for persons in motion.

3 MSE

 3 wavelettype ( p) (q )

Fig. 2 :

 2 Fig. 2: Front side measurement: wavelet decomposition signal at level 8 (d8) using Coif3 Vs the corresponding ECG signal Fig. 3 shows the heartbeat signal extracted from the cardiopulmonary signal measured from the back side.

Fig. 3 :

 3 Fig. 3: Back side measurement: wavelet decomposition signal at level 8 (d8) using Coif3 Vs the corresponding ECG signal

Fig. 4 :

 4 Fig. 4: Left side measurement: wavelet decomposition signal at level 8 (d8) using Coif3 Vs the corresponding ECG signal

Fig. 5 :

 5 Fig. 5: Right side measurement: wavelet decomposition signal at level 8 (d8) using Coif3 Vs the corresponding ECG signal

Table 1 .

 1 HR absolute relative error using different wavelet families and for cardiopulmonary signals measured from the different sides of the subject

	Side	Frequency (GHz)	Bior 2.4	Rbior 1.3	Wavelet Family Sym5 Db5	Coif3	Dmey
		2.4	12	14	3	3	6	2
	Front	5.8	9	2	7	5	15	18
		10	8	4	7	4	3	17
		2.4	2	11	11	11	8	15
	Back	5.8	6	10	23	10	6	14
		10	13	12	8	17	15	17
		2.4	5	1	1	4	3	6
	Left	5.8	2	3	3	5	3	6
		10	8	2	2	3	9	6
		2.4	6	1	9	1	11	13
	Right	5.8	10	13	15	18	10	13
		10	15	17	15	15	13	9