A Gaudillière
email: gaudilli@cmi.univ-mrs.fr

J Reygner
email: julien.reygner@polytechnique.edu

Sampling the Fermi statistics and other conditional product measures 1

Keywords: Metropolis algorithm, Markov chain, sampling, mixing time, product measure, conservative dynamics AMS 2010 classification: 60J10, 82C44

Through a Metropolis-like algorithm with single step computational cost of order one, we build a Markov chain that relaxes to the canonical Fermi statistics for k non-interacting particles among m energy levels. Uniformly over the temperature as well as the energy values and degeneracies of the energy levels we give an explicit upper bound with leading term km ln k for the mixing time of the dynamics. We obtain such construction and upper bound as a special case of a general result on (nonhomogeneous) products of ultra log-concave measures (like binomial or Poisson laws) with a global constraint. As a consequence of this general result we also obtain a disorder-independent upper bound on the mixing time of a simple exclusion process on the complete graph with site disorder. This general result is based on an elementary coupling argument, illustrated in a simulation appendix and extended to (non-homogeneous) products of log-concave measures.

1 From the Fermi statistics to general conditional products of log-concave measures the canonical Fermi statistics at inverse temperature β for k non-interacting particles among the m energy levels 1, . . . , m, with energy values v 1 , . . . , v m and degeneracies n 1 , . . . , n m is the conditional probability measure on

X k,m := {(k 1 , . . . , k m ) ∈ N m : k 1 + • • • + k m = k} (1.2)
given by

ν := µ(•|X k,m ) (1.3)
with µ the product measure on N m such that

µ(k 1 , . . . , k m ) := 1 Z m j=1 n j k j exp{-βk j v j }, (1.4)
Z := k1,...,km m j=1 n j k j exp{-βk j v j },

with nj kj = 0 whenever k j > n j . In other words, ν is a (non-homogeneous) product of binomial laws in k 1 , . . . , k m with the global constraint

k 1 + • • • + k m = k (1.6)
and we can write

ν(k 1 , . . . , k m ) = 1 Q m j=1 e φj (kj ) , (k 1 , . . . , k m ) ∈ X k,m (1.7)
where the φ j are defined by

φ j : k j ∈ N → -βk j v j + ln n j k j ∈ R ∪ {-∞} (1.8)
and Q is such that ν is a probability measure.

The first aim of this paper is to describe an algorithm that simulates a sampling according to ν in a time that can be bounded from above by an explicit polynomial in k and m, uniformly over β, (v j ) 1≤j≤m and (n j ) 1≤j≤m .

The reason why we prefer a bound in k and m rather than in the 'volume' of the system n = j n j , will be clarified later.

A first naive (and wrong) idea to do so consists in choosing the position (the energy level) of a first, second, . . . and eventually k th particle in the following way. First choose randomly the position of the first particle according to the exponential weights associated with the 'free entropies' of the empty sites, that is choose level j with a probability proportional to exp{-βv j + ln n j }. Then decrease by 1 the degeneracy of the chosen energy level and repeat the procedure to choose the position of the second, third, . . . and eventually k th particle. It is easy to check that, doing so, the final distribution of the occupation numbers k 1 , . . . , k m associated with the different energy levels, that is of the numbers of particles placed in each level, is in general not given by ν as soon as k is larger than one. But it turns out that this naive idea can be adapted to build an efficient algorithm to perform approximate samplings under the Fermi statistics.

Very classically, the fast sampling performed by the algorithm we will build will be obtained by running a Markov chain X with transition matrix p on X k,m and with equilibrium measure ν. The efficiency of the algorithm will be measured through the bounds that we will be able to give on the mixing time t ǫ , defined for any positive ǫ < 1 by As a consequence, estimating mixing times is not the only one issue of this paper, building a 'good' Markov chain is part of the problem. As far as that part of the problem is concerned, we propose to build a Metropolis-like algorithm that uses the 'free energies' of the naive approach to define a conservative dynamics. Assuming that at time t ∈ N the system is in some configuration X t = η in X k,m with ν(η) > 0, and defining for any η = (k 1 , . . . , k m ) and any distinct i and j in {1; • • • ; m}

t ǫ := inf{t ≥ 0 : d(t) ≤ ǫ} (1.9) d(t) := max
η ij := (k ′ 1 , . . . , k ′ m ) with k ′ s =    k s for s ∈ {1; • • • ; m} \ {i; j} k i -1 if s = i k j + 1 if s = j ,
(1.12) the configuration at time t + 1 will be decided as follows:

• choose a particle with uniform probability (it will stand in a given level i with probability k i /k),

• choose an energy level with uniform probability (a given level j will be chosen with probability 1/m),

• with i the level where stood the chosen particle and j the chosen energy level, extract a uniform variable U on [0; 1) and set X t+1 = η ij if i = j and U < exp{-βv j + ln(n j -k j ) + βv iln(n i -(k i -1))}, (1.13)

X t+1 = η if not.
In other words, denoting by [a] + = (a + |a|)/2 the positive part of any real number a and with

ψ j : k j ∈ {0; • • • ; n j } → -βv j + ln(n j -k j ) ∈ R ∪ {-∞}, j ∈ {1; • • • ; m}, (1.14) for any distinct i and j in {1; • • • ; m} P (X t+1 = η ij |X t = η) = p(η, η ij ) = k i k 1 m exp{-[ψ i (k i -1) -ψ j (k j )] + }, (1.15) and P (X t+1 = η|X t = η) = p(η, η) = 1 - i =j p(η, η ij ). (1.16)
Remark: In order to avoid any ambiguity in (1.15) in the case k i = 0, we set ψ i (-1) = +∞ (even though the algorithm we described does not require any convention for ψ i (-1)). This Markov chain is certainly irreducible and aperiodic. To prove that it relaxes to ν we will check the reversibility of the process with respect to ν. Then we will have to estimate the mixing time of the process. We will carry out both the tasks in a more general setup.

A general result

For any function f : N → R we define

∇ + f : x ∈ N → ∇ + x f := f (x + 1) -f (x) (1.17) ∇ -f : x ∈ N \ {0} → ∇ - x f := f (x -1) -f (x) (1.18) ∆f : x ∈ N \ {0} → ∆ x f := ∇ + x f + ∇ - x f = -∇ - x (∇ + f ) (1.19)
and we say that a measure γ on the integers

γ : x ∈ N → e φ(x) ∈ R + , (1.20) with φ : N → R ∪ {-∞}, is log-concave if N \ γ -1 ({0}) is an interval of the integers and γ(x) 2 ≥ γ(x -1)γ(x + 1), x ∈ N \ {0}, (1.21)
i.e., if ∇ + φ is non-increasing, or, equivalently, -∆φ is non-negative (with the obvious extension of the previous definitions to such a possibly non-finite φ).

The measure µ defined in (1.4) is a product of log-concave measures and the canonical Fermi statistics is such a product measure normalized over the condition (1.6).

N.B.

From now on, and except for explicit mentioning of additional hypotheses, we will only assume that the probability measure ν we want to sample is a product of log-concave measures normalized over the global constraint (1.6), i.e., that ν is a probability on X k,m that can be written in the form (1.7) with non-increasing ∇ + φ j 's.

In this more general setup we will often refer to the indices j in {1; • • • ; m} as sites rather than energy levels of the system.

Actually the e φj 's of the Fermi statistics are much more than log-concave measures. They are ultra log-concave measures according to the following definition by Pemantle [START_REF] Pemantle | Towards a theory of negative dependence[END_REF] and Liggett [START_REF] Liggett | Ultra logconcave sequences and negative dependence[END_REF].

Definition 1.2.1 A measure γ : N → R + is ultra log-concave if x → x!γ(x) is log-concave.
In other words e φj is ultra log-concave if and only if

ψ j := ∇ + φ j + ln(1 + •) (1.22)
is non-increasing (for the Fermi statistics observe that so are the φ j 's and that (1.22) is consistent with (1.14)).

For birth and death processes that are reversible with respect to ultra logconcave measures, Caputo, Dai Pra and Posta [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF] proved modified log-Sobolev inequalities and stronger convex entropy decays, both giving good upper bounds on the mixing time of the processes. Johnson [START_REF] Johnson | Bounds on the Poincaré constant of ultra log-concave random variables[END_REF] proved also easier Poincaré inequalities that give weaker bounds on the mixing times. We refer to [START_REF] Saloff-Coste | Lectures on finite Markov chains[END_REF], [START_REF] Ané | Sur les inégalités de Sobolev logarithmiques Panoramas et Synthèses no[END_REF], [START_REF] Montenegro | Mathematical Aspects of Mixing Times in Markov Chains[END_REF] for an introduction to this classical functional inequality approach to convergence to equilibrium and we note that for such birth and death processes the ultra-log concavity hypotheses allowed for a Bakry-Émery like approach (see [START_REF] Bakry | Diffusions hypercontractives, Séminaire de probabilités[END_REF]) to derive (modified) log-Sobolev inequalities. Actually, [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF] was an attempt to extend this celebrated analysis to Markov processes with jumps (see also [START_REF] Joulin | Curvature, concentration, and error estimates for Markov chain Monte Carlo[END_REF] for a more geometric perspective). But it turns out, as we will soon review, that beyond the case of birth and death processes the role of ultra logconcavity to follow the Bakry-Émery line in the discrete setup is still unclear.

In this paper we will not follow the functional analysis approach to control mixing times. To bound from above the mixing time of a Markov chain that is reversible with respect to a conditional product of ultra log-concave measures, we will follow (and recall in Section 2) the non less classical, and, in this case, elementary, probabilistic approach via coalescent coupling. We will prove:

Theorem 1 If ν derives from a product of ultra log-concave measures, then the Markov chain with transition matrix p defined by

p(η, η ij ) = k i k 1 m exp{-[ψ i (k i -1) -ψ j (k j )] + }, η = (k 1 , . . . , k m ) ∈ X k,m \ ν -1 {0}, i = j ∈ {1; • • • ; m}, (1.23) p(η, η) = 1 - i =j p(η, η ij ) (1.24)
is reversible with respect to ν and, for any positive ǫ < 1, its mixing time t ǫ satisfies t ǫ ≤ km ln(k/ǫ).

(1.25)

Proof: see Section 2.
The most relevant point of Theorem 1 with respect to the previous results we know stands in the uniformity of the upper bound above the disorder of the system (except for the ultra log-concavity hypothesis on e φj in each j). In particular and as far as the Fermi statistics is concerned, our estimate does not depend on the temperature, and, more generally it is independent from the energy values as well as the level degeneracies.

To illustrate this fact let us start with the case n j = 1 for all j. In this case our dynamics is a simple exclusion process with site disorder. Caputo ([6], [START_REF] Caputo | On the spectral gap of the Kac walk and other binary collision processes[END_REF]) proved Poincaré inequalities for such processes, in their continuous time version, assuming a uniform lower (and upper) bound on general transition rates, while Caputo, Dai Pra and Posta [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF], looking at particular rates for the process and still assuming moderate disorder -that is, uniform lower and upper bounds on these rates -proved a modified log-Sobolev inequality. For the particular choice of rates they made, the upper bound on the mixing time implied by [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF] could not hold in a strong disorder context (for example with k = 1, m = 3, v 1 = v 2 = 0, v 3 > 0 and β ≫ 1). Our uniformity over the disorder of the system depends then strongly on our particular choice for the transition probabilities. As it is often the case with Markov processes on discrete state space, the details of the dynamics are not less important than the properties of its equilibrium measure.

Under the moderate disorder hypotheses of [START_REF] Caputo | Spectral gap inequalities in product spaces with conservation laws[END_REF], [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF], [START_REF] Caputo | On the spectral gap of the Kac walk and other binary collision processes[END_REF], however, the upper bound on the mixing time implied by [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF] and suggested by [START_REF] Caputo | Spectral gap inequalities in product spaces with conservation laws[END_REF], [START_REF] Caputo | On the spectral gap of the Kac walk and other binary collision processes[END_REF] is better than our upper bound in Theorem 1 (by a factor of order k). But we will see that introducing such moderate disorder hypotheses in our arguments directly improves our result by a factor of order m ≥ k (see the last remark at the end of Section 2).

To close the discussion on simple exclusion processes we note that those of the arguments in [START_REF] Caputo | On the spectral gap of the Kac walk and other binary collision processes[END_REF] that do not depend on the disorder suggest an upper bound on the mixing time of order m 2 ln m. Theorem 1 improves this estimate when k is small with respect to m.

As far as more general conditional non-homogeneous product of log-concave measures are concerned, we stress once again that Theorem 1 gives a uniform bound over the disorder that can be improved by a factor of order m by adding moderate disorder hypotheses (see the last remark at the end of Section 1.4 and our simulation Appendix). Then, such product measures can be equilibrium measures of zero-range processes with a continuous time generator as in [START_REF] Boudou | Spectral gap inequalities for interacting particle systems via a Bochner type inequality[END_REF], [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF]:

L η f := 1 m i =j c i (k i ) f (η ij ) -f (η) . (1.26)
where the c j : k j ∈ N → [0, +∞) are such that c j (0) = 0 and c j (k j ) > 0 for k j > 0. Indeed, such a process is reversible with respect to ν provided that

∀j ∈ {1; • • • ; m}, ∀k j ∈ N, φ j (k j ) = 0<l≤kj ln 1 c j (l)
.

(1.27) Boudou, Caputo, Dai Pra and Posta [START_REF] Boudou | Spectral gap inequalities for interacting particle systems via a Bochner type inequality[END_REF] proved a Poincaré inequality for such a process, assuming that there exists a positive c such that

∀j ∈ {1; • • • ; m}, ∀k j ∈ N, c j (k j + 1) -c j (k j ) ≥ c. (1.28)
This is a moderate disorder hypothesis that implies the log-concavity of the µ j .

To go to modified log-Sobolev inequalities, i.e., to good mixing time estimates rather than simple gap estimates, there is the additional hypothesis in [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF] that there exists a non-negative δ < c such that ∀j ∈ {1; • • • ; m}, ∀k j ∈ N, c j (k j + 1) -c j (k j ) ≤ c + δ.

(1.29)

Clearly there are ultra log-concave measures that do not satisfy (1.29): to have an ultra log-concave measure one needs a strongly decreasing ∇ + φ j , i.e., a strongly increasing c j . Conversely, it is not true that (1.28) and (1.29) imply ultra log-concavity, i.e.,

∀j ∈ {1; • • • ; m}, ∀k j ≥ 1, ln k j + 1 c j (k j + 1) ≤ ln k j c j (k j )
.

(1.30) However, elementary algebra shows that (1.28) together with (1.29) implies

∀j ∈ {1; • • • ; m}, ∀k j ≥ 1, ln k j + 1/2 c j (k j + 1) ≤ ln k j c j (k j ) (1.31)
and this strangely looks like (1.30). Therefore we said that the role played by ultra log-concavity is still unclear along the functional analysis line of research in the discrete setup. We conclude stressing once again that in this paper we will follow a different line, that our uniform bound on the mixing time comes from an elementary coupling argument and that we do not need any moderate disorder hypothesis.

Interpolating between sites and particles

It seems that today available techniques are such that the less log-concavity we have, the more homogeneity we need to control the convergence to equilibrium. Staying to the papers mentioned above, without ultra log-concavity or at least something that looks like ultra log-concavity we only have Poincaré inequalities for non-homogeneous product of log-concave measures, and without log-concavity we have modified log-Sobolev inequalities for homogeneous product measures only (see [START_REF] Caputo | Convex entropy decay via the Bochner-Bakry-Emery approach[END_REF]). In addition, the only result we know for a conservative dynamics in (weakly) disordered context and with an equilibrium measure that is a product of measures that are not log-concave is that of Landim and Noronha Neto [START_REF] Landim | Poincaré and logarithmic Sobolev inequality for Ginzburg-Landau processes in random environment[END_REF] for the (continuous) Ginzburg-Landau process.

We will see that all the ideas of the proof of Theorem 1 can be extended to deal with a large class of conditional product of log-concave measures that are not ultra log-concave. To do so, let us define

δ := max λ ∈ [0; 1] : ∀j ∈ {1; • • • ; m}, ∀k j > 0, -∆ kj φ j ≥ λ ln 1 + k j k j .
(1.32) In other words, δ is the largest real number in [0; 1] for which all the

ψ [δ] j := ∇ + φ j + δ ln(1 + •), j ∈ {1; • • • ; m}, (1.33)
are non-increasing. Denoting by a ∧ b the minimum of two real numbers a, b and defining

l δ := k δ (k ∧ m) 1-δ (1.34)
we will prove:

Theorem 2 If δ > 0 then the Markov chain with transition matrix p defined by

p(η, η ij ) = k δ i l δ 1 m exp -ψ [δ] i (k i -1) -ψ [δ] j (k j ) + , η = (k 1 , . . . , k m ) ∈ X k,m \ ν -1 {0}, i = j ∈ {1; • • • ; m}, (1.35) p(η, η) = 1 - i =j p(η, η ij ) (1.36)
is reversible with respect to ν and, for any positive ǫ < 1, its mixing time t ǫ satisfies

t ǫ ≤ (k ∧ m) 1-δ δ km ln(k/ǫ). (1.37) Remark 1: By Hölder's inequality, if δ > 0 then m i=1 k δ i = m i=1 k δ i 1 l 1-δ {ki =0} ≤ m i=1 k i δ m i=1 1 l {ki =0} 1-δ ≤ l δ (1.38)
and this ensures that (1.35)-(1.36) define a probability matrix. Remark 2: As far as this can make sense in our discrete setup, we note that the hypothesis δ > 0 is slightly weaker than a "uniform strict log-concavity hypothesis" (see (1.32)).

Proof of Theorem 2: see Section 3.

For δ = 1 the transition matrix represents an algorithm starting with a uniform choice of a particle. For δ = 0 Theorem 2 is empty but (1.35)-(1.36) still define a Markov chain X that can be seen as a particular version of a discrete state space non-homogeneous Ginzburg-Landau process. In this case the transition matrix represents an algorithm that starts with a uniform choice of a (non-empty) site. The case 0 < δ < 1 can be seen as an interpolation between uniform choices of site and particle. More precisely, assuming that at time t ∈ N the system is in some configuration

X t = η = (k 1 , . . . , k m ) in X k,m \ ν -1 {0}
, the configuration at time t + 1 can be decided as follows.

• Choose a site i or no site at all with probabilities proportional to k δ i and l δi k δ i . • If some site i was chosen, then proceed as in the previously described algorithm using the functions ψ

[δ] j instead of the ψ j 's, if not, then set X t+1 = η.

Last remarks and original motivation

First, we note that as long as one wants bounds that are uniform over the disorder, Theorem 1 gives the right order for the mixing time: for the Fermi statistic in the very low temperature regime with, for example, v 1 < v 2 < • • • < v m and n 1 , n 2 ≥ k, the equilibrium measure will be concentrated on (k, 0, . . . , 0) while, starting from (0, k, 0, . . . , 0), the system will reach the ground state in a time of order km ln k for large k (this is a coupon-collector estimate). This fact is illustrated in our simulation Appendix.

Next, we observe that, with our definitions, p(η, η) can often be close to one, especially in strong disorder situations or when the right and left hand sides in (1.38) are far from each other. If one would like to use these results to perform practical simulations, then it could be useful to note that the computational time would still be improved by implementing an algorithm that at each step simulates, for a given configuration η on the trajectory of the Markov chain, the elapsed time before the particle reach a different configuration (this is a geometric time) and choose this configuration η ′ = η according to the (easy to compute) associated law. It would then be enough to stop the algorithm as soon as the total simulated time goes beyond the mixing time (and then return the last configuration, that the system occupied at the mixing time) rather than waiting for the original algorithm to make a step number equal to the mixing time.

Turning back to the first naive and wrong idea, it is interesting to note that it can easily be modified to determine the most probable states of the system, i.e., the configurations η

* = (k * 1 , . . . , k * m ) in X k,m such that m j=1 φ j (k * j ) = max k1+•••+km=k m j=1 φ j (k j ). (1.39)
One can prove, using the concavity of the φ j 's, that the most probable configurations for the system with k particles can be obtained from the most probable configurations (k ′ 1 , . . . , k ′ m ) for the system with k -1 particles simply adding one particle where the corresponding gain in 'free energy' is the highest, that is in j * such that

∇ + k ′ j * φ j * = max j ∇ + k ′ j φ j . (1.40)
As a consequence one can place the particles one by one, each time maximizing this free energy gain, to build the most probable configuration. Then, as a referee pointed out, since the sampling problem is a trivial one in the Poissonian case (when all the µ j 's are Poisson measures, so that ν is nothing but a multinomial law M ), one can ask about the expected time needed to perform a rejection sampling with respect to the multinomial case. It is given by max η ν(η)/M (η). If we take the Fermi statistics with m = k and all the n j equal to 1, then we find (optimizing on the multinomial parameters by a geometric/arithmetic mean comparison)

k k /k! ∼ e k / √ 2πk.
Of course, the sampling problem for that Fermi statistics is also a trivial one. But if we take β = 0 and all the n j equal to 2, it is not anymore a trivial problem and we find an expected time for the rejection sampling that is logarithmically equivalent to (e/2) k .

Turning back to the Fermi statistics we now explain why we were interested in bounds in k and m rather than in the volume n. Iovanella, Scoppola and Scoppola defined in [START_REF] Iovanella | Some spin glass ideas applied to the clique problem[END_REF] an algorithm to individuate cliques (i.e., complete subgraphs) with k vertices inside a large Erdös-Reyni random graph with n vertices. Their algorithm requires to perform repeated approximate samplings of Fermi statistics in volume n, with k particles and m = 2k + 1 energy levels. Now, the key observation is that the largest cliques in Erdös-Reyni graphs with n vertices are of order ln n, so that k and m in this problem are logarithmically small with respect to n. Before Theorem 1 the samplings for their algorithm were done by running simple exclusion processes with k particles on the complete graph (with site disorder) of size n. Such processes converge to equilibrium in a time of order kn ln k. Now the samplings are done in a time of order km ln k ∼ 2k 2 ln k, and that was the original motivation of the present work.

Finally we note that our bound in k and m is not only useful when k is small with respect to n but also when k is close to n. In this case we can define a dynamics on the n -k vacancies rather than on the k particles.

Proof of Theorem 1

In this section we assume that ν is a measure on X k,m deriving from a product of ultra log-concave measures, which means that we can write ν(k 1 , . . . , k m ) = (1/Q) j e φj (kj ) and, for all j ∈ {1; • • • ; m}, ψ j defined by (1.22) is nonincreasing.

Reversibility

We first prove that the transition matrix p defined by (1.23) and (1.24) is reversible with respect to the measure ν.

Let η = (k 1 , . . . , k m ) ∈ X k,m \ ν -1 {0} and let i = j ∈ {1; • • • ; m}. We have p(η, η ij ) = 0 if and only if ν(η ij ) = 0 and, in that case, ν(η ij ) ν(η) = e φi(ki-1)+φj (kj +1) e φi(ki)+φj (kj ) = exp{∇ - ki φ i + ∇ + kj φ j } (2.1) while p(η ij , η) p(η, η ij ) = k j + 1 k i exp{-[ψ j (k j ) -ψ i (k i -1)] + +[ψ i (k i -1) -ψ j (k j )] + } (2.2) = exp{ψ i (k i -1) -ψ j (k j ) + ln(k j + 1) -ln(k i )} (2.3) = exp{∇ + ki-1 φ i -∇ + kj φ j } (2.4) = exp{-∇ - ki φ i -∇ + kj φ j } (2.5) so that ν(η)p(η, η ij ) = ν(η ij )p(η ij , η). (2.6)

A few words about the coupling method

In order to upper bound the mixing time of the Markov chain with transition matrix p, we will use the coupling method. Given a Markov chain (X 1 , X 2 ) on X k,m × X k,m , we say it is a (Markovian) coupling for the dynamics if both X 1 and X 2 are Markov chains with transition matrix p. Given such a coupling, we define the coupling time τ couple as the first (random) time for which the chains meet, that is

τ couple := inf{t ≥ 0 : X 1 t = X 2 t }. (2.7)
In this work, every coupling will also satisfy the condition

t ≥ τ couple ⇒ X 1 t = X 2 t . (2.8)
Then, it is a well-known fact that for all t ≥ 0,

d(t) ≤ max η,θ∈X k,m P (τ couple > t|X 1 0 = η, X 2 0 = θ) (2.9)
where d(t) is defined by (1.10). A proof of this fact as well as an exhaustive introduction to mixing time theory can be found in [START_REF] Levin | Markov Chains and mixing times[END_REF].

In the proof of both 1 and 2 we will build a coupling for which there exists a function ρ that measures in some sense a 'distance' between X 1 t and X 2 t and from which we will get a bound on the mixing time thanks to the following proposition.

Proposition 2.2.1 Let (X 1 , X 2 ) be a coupling for a Markov chain with transition matrix p. We assume that the coupling satisfies the property (2.8). Let ρ : X k,m × X k,m → N such that ρ(η, θ) = 0 if and only if η = θ. If M is the maximum of ρ and if there exists α > 1 such that, for all t ≥ 0,

E ρ(X 1 t+1 , X 2 t+1 )|X 1 t , X 2 t ≤ 1 - 1 α ρ(X 1 t , X 2 t ), (2.10)
then for all ǫ > 0 the mixing time t ǫ of the dynamics is upper bounded by α ln(M/ǫ).

Proof: Remark that (2.10) actually means that (ρ(X

1 t , X 2 t )(1 -1/α) -t
) t∈N is a super-martingale. Taking the expectation we get

E(ρ(X 1 t+1 , X 2 t+1 )) ≤ 1 - 1 α E(ρ(X 1 t , X 2 t )) (2.11)
As a consequence

E(ρ(X 1 t , X 2 t )) ≤ 1 - 1 α t E(ρ(X 1 0 , X 2 0 )) ≤ M e -t α .
(2.12)

Since we assume ρ(η, θ) = 0 if and only if η = θ and (2.8),

P (τ couple > t) = P (ρ(X 1 t , X 2 t ) > 0) = P (ρ(X 1 t , X 2 t ) ≥ 1). (2.13)
From Markov's inequality and (2.12) we deduce

P (τ couple > t) ≤ E(ρ(X 1 t , X 2 t )) ≤ M e -t α .
(2.14)

Since the upper bound in (2.14) is uniform in X 1 0 and X 2 0 , according to (2.9) we get d(t) ≤ M e -t α .

(2.15)

Thus, given ǫ > 0, if t ≥ α ln(M/ǫ) then d(t) ≤ ǫ, so that t ǫ ≤ α ln(M/ǫ).

A colored coupling

In (2.17)

We will couple two dynamics ω 1 and ω 2 on Ω, then we will work with the coupling (ξ(ω 1 ), ξ(ω 2 )). We will build the coupling (ω 1 , ω 2 ) thanks to the coloring we now introduce. At step t, a red-blue coloring of (ω 1 t , ω 2 t ) is a couple of functions C 1 , C 2 : {1; . . . ; k} → {blue; red} such that for all energy level i ∈ {1; . . . ; m}, the number of blue particles in the level i is the same in both distributions, i.e., if |A| refers to the cardinality of the set A,

|{x : C 1 (x) = blue, ω 1 t (x) = i}| = |{x : C 2 (x) = blue, ω 2 t (x) = i}| (2.18)
and an energy level cannot contain red particles in both distributions

C 1 (x) = red ⇒ ∀y, (ω 2 t (y) = ω 1 t (x) ⇒ C 2 (y) = blue), (2.19) C 2 (x) = red ⇒ ∀y, (ω 1 t (y) = ω 2 t (x) ⇒ C 1 (y) = blue). (2.20)
As a consequence of (2.18), there exists a one-to-one correspondence

Φ : {x : C 1 (x) = blue} → {x : C 2 (x) = blue} (2.21)
such that for all x, ω 2 t (Φ(x)) = ω 1 t (x). Since the number of blue particles is the same in both distributions, so is the number of red particles. Moreover, ξ(ω 1 t ) = ξ(ω 2 t ) if and only if all particles are blue. Then we are willing to provide a coupling (ω 1 , ω 2 ) for which the number of red particles is non-increasing. This number does not depend on the red-blue coloring and at any time t we will call it ρ t . Writing ξ(ω

1 t ) = (k 1 1 , . . . , k 1 m ) and ξ(ω 2 t ) = (k 2 1 , . . . , k 2 m )
we have the identities

ρ t = 1 2 m i=1 |k 1 i -k 2 i | = m i=1 [k 1 i -k 2 i ] + = m i=1 [k 2 i -k 1 i ] + . ( 2

.22)

We are now ready to build our coupling (ω 1 , ω 2 ). Given (ω 1 t , ω 2 t ) and a redblue coloring (C 1 , C 2 ) at step t (such a coloring certainly exists for any couple of distributions (ω 1 t , ω 2 t )), let x 1 be a uniform random integer in {1; . . . ; k}.

• If C 1 (x 1 ) = blue: then we set x 2 = Φ(x 1 ) where Φ is provided by (2.21).

• If C 1 (x 1 ) = red: let x 2 be a uniform random integer in {x : C 2 (x) = red}.

Lemma 2.3.1 The random integer x 2 has a uniform distribution on the set {1; . . . ; k}.

Proof: We write

P (x 2 = y) = k x=1 P (x 2 = y|x 1 = x)P (x 1 = x) (2.23) = x:C 1 (x)=blue 1 l {y=Φ(x)} 1 k + x:C 1 (x)=red 1 l {C 2 (y)=red} ρ t 1 k (2.24) = 1 l {C 2 (y)=blue} k + 1 l {C 2 (y)=red} ρ t ρ t k (2.25) = 1 k . ( 2

.26)

We then choose an energy level j uniformly in {1; . . . ; m} and we write ξ(ω

1 t ) = (k 1 1 , . . . , k 1 m ) and ξ(ω 2 t ) = (k 2 1 , . . . , k 2 m ).
• If C 1 (x 1 ) = blue: then C 2 (x 2 ) = blue, and x 1 , x 2 are in the same energy level

i := ω 1 t (x 1 ) = ω 2 t (x 2 ). Let a = b ∈ {1; 2} such that ψ i (k a i -1) - ψ j (k a j ) ≤ ψ i (k b i -1) -ψ j (k b j ). Then, p a := exp{-[ψ i (k a i -1) -ψ j (k a j )] + } (2.27) ≥ exp{-[ψ i (k b i -1) -ψ j (k b j )] + } =: p b . (2.28)
Let U be a uniform random variable on [0; 1).

-If U < p b : then we set ω a t+1 (x) = ω a t (x) for all x = x a , ω a t+1 (x a ) = j, ω b t+1 (x) = ω b t (x) for all x = x b and ω b t+1 (x b ) = j. Then ξ(ω a t+1 ) = (ξ(ω a t )) ij and ξ(ω b t+1 ) = (ξ(ω b t )) ij , and for any red-blue coloring of (ω 1 t+1 , ω 2 t+1 ) the number of red particles ρ t+1 remains the same (since both particles x 1 and x 2 have moved together).

-If p b ≤ U < p a : then we set ω a t+1 (x) = ω a t (x) for all x = x a , ω a t+1 (x a ) = j and ω b t+1 (x) = ω b t (x) for all x. Then ξ(ω a t+1 ) = (ξ(ω a t )) ij and ξ(ω b t+1 ) = ξ(ω b t ). The only situation in which the number of red particles could increase is the following: k a j ≥ k b j and k a i ≤ k b i . Since ψ i and ψ j are non-increasing, this would imply p a ≤ p b that contradicts p b ≤ U < p a . Then, ρ t+1 ≤ ρ t for any red-blue coloring of (ω

1 t+1 , ω 2 t+1 ). -If U ≥ p a : then we set ω a t+1 = ω a t and ω b t+1 = ω b t . Then ξ(ω a t+1 ) = ξ(ω a t ), ξ(ω b t+1 ) = ξ(ω b t ) and ρ t+1 = ρ t . • If C 1 (x 1 ) = red: then C 2 (x 2 ) =
red, and we define i 1 := ω 1 t (x 1 ) and i 2 := ω 2 t (x 2 ). Note that according to (2.19), i 1 = i 2 . Let then V 1 , V 2 be two independent uniform random variables on [0; 1).

-If V 1 < exp{-[ψ i 1 (k 1 i 1 -1) -ψ j (k 1 j )] + } then we set ω 1 t+1 (x) = ω 1 t (x) for all x = x 1 and ω 1 t+1 (x 1 ) = j and then ξ(ω 1 t+1 ) = (ξ(ω 1 t )) ij ; otherwise we leave ω 1 t+1 = ω 1 t . -If V 2 < exp{-[ψ i 2 (k 2 i 2 -1) -ψ j (k 2 j )
] + } then we set ω 2 t+1 (x) = ω 2 t (x) for all x = x 2 and ω 2 t+1 (x 2 ) = j and then ξ(ω 2 t+1 ) = (ξ(ω 2 t )) ij ; otherwise we leave ω 2 t+1 = ω 2 t .

Whether red particles move or not, the number of red particles cannot increase, so it is clear that ρ t+1 ≤ ρ t .

We conclude Proposition 2.3.2 (ρ t ) t∈N is a non-increasing process.

and claim

Proposition 2.3.3 The process (ξ(ω 1 ), ξ(ω 2 )) is a coupling for the Markov chain with transition matrix p.

Proof: Writing ξ(ω 1 t ) = η and given i, j ∈ {1; . . . ; m},

P (ξ(ω 1 t+1 ) = η ij ) = P (ξ(ω 1 t+1 ) = η ij |C 1 (x 1 ) = blue)P (C 1 (x 1 ) = blue) + P (ξ(ω 1 t+1 ) = η ij |C 1 (x 1 ) = red)P (C 1 (x 1 ) = red). (2.29)
We have

P (C 1 (x 1 ) = red) = 1 -P (C 1 (x 1 ) = blue) = ρ t k (2.30)
and

P (ξ(ω 1 t+1 ) = η ij |C 1 (x 1 ) = blue) = P (ω 1 t (x 1 ) = i|C 1 (x 1 ) = blue) × 1 m × P (U < p 1 ) (2.31) = k 1 i ∧ k 2 i (k -ρ t )m exp{-[ψ i (k 1 i -1) -ψ j (k 1 j )] + } (2.32) P (ξ(ω 1 t+1 ) = η ij |C 1 (x 1 ) = red) = P (ω 1 t (x 1 ) = i|C 1 (x 1 ) = red) × 1 m × P (V 1 < p 1 ) (2.33) = [k 1 i -k 2 i ] + ρ t m exp{-[ψ i (k 1 i -1) -ψ j (k 1 j )] + } (2.34)
which finally leads to

P (ξ(ω 1 t+1 ) = η ij ) = k 1 i km exp{-[ψ i (k 1 i -1) -ψ j (k 1 j )] + } = p(η, η ij ). (2.35)
Then, P (ξ(ω 1 t+1 ) = η ij |ω 1 t , ω 2 t ) depends on ξ(ω 1 t ) only, which means that ξ(ω 1 ) is a Markov chain. Besides, according to (2.35) its transition matrix is p. Finally, by Lemma 2.3.1, the same is true for ξ(ω 2 ).

From now on, we will write X 1 t := ξ(ω 1 t ) and X 2 t := ξ(ω 2 t ). The previous proposition ensures that (X 1 , X 2 ) is a coupling for the dynamics with transition matrix p.

Estimating the coupling time

We will use Proposition 2.2.1. Since (ρ t ) t∈N is non increasing and ρ t = 0 if and only if X 1 t = X 2 t , all we have to do is to estimate from below the probability of the event {ρ t+1 < ρ t }.

Proposition 2.4.1 If at step t of the coupled dynamics (ω 1 , ω 2 ) we assume that red particles have been chosen, i.e., C 1 (x 1 ) = red = C 2 (x 2 ), then there is a choice of j ∈ {i 1 ; i 2 } for which the number of red particles decreases with probability 1, where i 1 (resp. i 2 ) still refers to ω 1 t (x 1 ) (resp. ω 2 t (x 2 )).

Proof:

If the inequalities exp{-[ψ i 1 (k 1 i 1 -1) -ψ i 2 (k 1 i 2 )] + } < 1 (2.36) exp{-[ψ i 2 (k 2 i 2 -1) -ψ i 1 (k 2 i 1 )] + } < 1 (2.37) hold together, then ψ i 1 (k 1 i 1 -1) > ψ i 2 (k 1 i 2 ) and ψ i 2 (k 2 i 2 -1) > ψ i 1 (k 2 i 1 ). Besides, since C 1 (x 1 ) = red, according to (2.19) k 1 i 1 > k 2 i 1 from which we get k 1 i 1 -1 ≥ k 2 i 1 and, since ψ i 1 is non-increasing, ψ i 1 (k 1 i 1 -1) ≤ ψ i 1 (k 2 i 1 ). Likewise, since C 2 (x 2 ) = red we have ψ i 2 (k 2 i 2 -1) ≤ ψ i 2 (k 1 i 2 )
. We finally may write

ψ i 1 (k 1 i 1 -1) > ψ i 2 (k 1 i 2 ) ≥ ψ i 2 (k 2 i 2 -1) > ψ i 1 (k 2 i 1 ) ≥ ψ i 1 (k 1 i 1 -1) (2.38)
which is absurd. As a result, either (2.36) or (2.37) is false. For instance let us assume that (2.36) is false, then if j = i 2 , with probability 1 we have

ω 1 t+1 (x) = ω 1 t (x) for all x = x 1 , ω 1 t+1 (x 1 ) = i 2 = ω 1 t (x 1
) and ω 2 t+1 = ω 2 t . Then, the number of red particles for any red-blue coloring of (ω 1 t+1 , ω 2 t+1 ) is exactly ρ t+1 = ρ t -1.

It follows

Corollary 2.4.2 At step t, if ρ t > 0, the probability for ρ t+1 to be ρ t -1 is at least ρ t /km. Consequently, and owing to the fact that ρ t cannot increase, we have the inequality

E(ρ t+1 -ρ t |X 1 t , X 2 t ) ≤ -P (ρ t+1 = ρ t -1|X 1 t , X 2
t ) from which we deduce

E(ρ t+1 |X 1 t , X 2 t ) ≤ 1 - 1 km ρ t .
(2.39)

Then we can apply Proposition 2.2.1 to ρ(X 1 t , X 2 t ) = ρ t with M = k and α = km, which finally proves Theorem 1.

Remark: Adding "moderate disorder hypotheses" we can gain a lot, depending on the specific model we consider. For example in the case of our simple exclusion dynamics, i.e., for the Markov chain associated with the Fermi statistics when all the n j 's are equal to 1, we can gain a factor of order m. We can indeed assume that k/m ≤ 1/2 (if not, we consider a dynamics on the vacancies rather than on the particles as mentioned in Section 1.4) and, if, for example, ν derives from an homogeneous product measure, i.e., ψ i = ψ j for all i and j, then there are much more than one choice for j for which the number of red particles will decrease with probability 1 when red particles were chosen: any of the more than m/2 vacant sites in X 1 t or X 2 t will do the job. In this case we gain a factor m/2, first in Corollary 2.4.2, then in the final estimate. More generally, if the functions |ψ i -ψ j | are uniformly bounded (this is a moderate disorder hypothesis), then, when red particles are chosen together with one of these more than m/2 vacant sites, the probability that the number of red particles decrease is bounded away from zero: once again we gain a factor of order m.

Proof of Theorem 2

We now work with the dynamics defined by (1.35)-(1.36) and we assume δ > 0. First, the reversibility of this dynamics with respect to ν still holds, with exactly the same computation as in Subsection 2.1. However, it is no longer possible to work with an underlying process ω t ∈ Ω since the factor k δ i cannot stand for a number of particles as soon as δ < 1. Therefore we need to adapt the coupling (X 1 , X 2 ) directly on X k,m .

Generalizing the previous coupling

At step t, let us assume

X 1 t = (k 1 1 , . . . , k 1 m ) ∈ X k,m and X 2 t = (k 2 1 , . . . , k 2 m ) ∈ X k,m .
We define the following sets:

R 1 := {i ∈ {1; . . . ; m} : k 1 i > k 2 i }, (3.1) R 2 := {i ∈ {1; . . . ; m} : k 2 i > k 1 i }, (3.2) B := {i ∈ {1; . . . ; m} : k 1 i = k 2 i }, (3.3)
and the following quantities:

w 1 := i∈R1 (k 1 i ) δ , w ′ 1 := i∈R2 (k 1 i ) δ , (3.4)
w 2 := i∈R2 (k 2 i ) δ , w ′ 2 := i∈R1 (k 2 i ) δ , (3.5)
w B := i∈B (k 1 i ) δ = i∈B (k 2 i ) δ . (3.6)
Finally we define

ρ t := 1 2 m i=1 |k 1 i -k 2 i | = m i=1 [k 1 i -k 2 i ] + = m i=1 [k 2 i -k 1 i ] + . (3.7)
Keeping in mind the previous coloring, R 1 (resp. R 2 ) is the set of sites in which there are red particles for the first (resp. the second) configuration, B is the set of sites in which there are only blue particles or no particles for both configurations, w 1 (resp. w 2 ) is proportional to the probability to choose a site for the first (resp. the second) configuration in which there are red particles, w ′ 1 (resp. w ′

2 ) is proportional to the probability to choose a site for the first (resp. the second) configuration in which there are only blue particles while there are red particles in the second (resp. the first) configuration, w B is proportional to the probability to choose a site in which there are only blue particles for both configurations, and l δ -

(w B + w 1 + w ′ 1 ) (resp. l δ -(w B + w 2 + w ′ 2 )
) is proportional to the probability not to choose any site for the first (resp. the second) configuration, and it can be positive as soon as δ < 1. Finally, ρ t still stands for the number of red particles and it is clear that ρ t = 0 if and only if

X 1 t = X 2 t . N.B.
In the remaining part of this subsection, we assume

w 1 + w ′ 1 ≥ w 2 + w ′ 2
in order to not overload the notations and not increase the number of cases to investigate. Obviously, the case w 1 + w ′ 1 ≤ w 2 + w ′ 2 is exactly symmetric. Let I be a uniform random variable on [0; l δ ).

(i) If I < w B : then there exists a unique i ∈ B such that

i ′ ∈B;i ′ <i (k 1 i ′ ) δ ≤ I < i ′ ∈B;i ′ ≤i (k 1 i ′ ) δ (3.8)
and we set i 1 = i 2 = i.

Remark: We then have, for all i ∈ {1; . . . ; m},

P (i 1 = i|(i)) = 1 l {i∈B} (k 1 i ) δ w B , (3.9)
P (i 2 = i|(i)) = 1 l {i∈B} (k 1 i ) δ w B = 1 l {i∈B} (k 2 i ) δ w B . (3.10) (ii) If w B ≤ I < w B + w ′ 1 : then there exists a unique i ∈ R 2 such that i ′ ∈R2;i ′ <i (k 1 i ′ ) δ ≤ I -w B < i ′ ∈R2;i ′ ≤i (k 1 i ′ ) δ (3.11)
and we set i 1 = i 2 = i.

Remark: We then have, for all i ∈ {1; . . . ; m},

P (i 1 = i|(ii)) = 1 l {i∈R2} (k 1 i ) δ w ′ 1 , (3.12)
P (i 2 = i|(ii)) = 1 l {i∈R2} (k 1 i ) δ w ′ 1 . (3.13) (iii) If w B + w ′ 1 ≤ I < w B + w ′ 1 + w 1 : then there exists a unique i ∈ R 1 such that i ′ ∈R1;i ′ <i (k 1 i ′ ) δ ≤ I -w B -w ′ 1 < i ′ ∈R1;i ′ ≤i (k 1 i ′ ) δ . (3.14)
We set i 1 = i and we define u

:= I -w B -w ′ 1 -i ′ ∈R1;i ′ <i (k 1 i ′ ) δ , so that 0 ≤ u < (k 1 i ) δ . Notice that since i ∈ R 1 , k 1 i > k 2 i . If u < (k 2 i ) δ then we set i 2 = i. Otherwise, for all i ′ ∈ R 2 we write k i ′ := ℓ∈R2;ℓ<i ′ (k 1
ℓ ) δ and we denote by T the disjoint union of intervals

T := i ′ ∈R2 k i ′ + (k 1 i ′ ) δ ; k i ′ + (k 2 i ′ ) δ ∪ [w 2 + w ′ 2 ; w 1 + w ′ 1 ). (3.15)
Let I ′ be a uniform random variable on T . If there exists

i ′ ∈ R 2 such that I ′ ∈ k i ′ + (k 1 i ′ ) δ ; k i ′ + (k 2 i ′ ) δ then we set i 2 = i ′ . Else we do not define i 2 .
Remark: We then have, for all i ∈ {1; . . . ; m},

P (i 1 = i|(iii)) = 1 l {i∈R1} (k 1 i ) δ w 1 (3.16)
and for all i ′ ∈ {1; . . . ; m},

P (i 2 = i ′ |(iii)) = 1 l {i ′ ∈R1} (k 2 i ′ ) δ w 1 + 1 l {i ′ ∈R2} i∈R1 (k 1 i ) δ -(k 2 i ) δ w 1 × (k 2 i ′ ) δ -(k 1 i ′ ) δ λ(T ) (3.17)
where

λ(T ) := i ′ ∈R2 (k 2 i ′ ) δ -(k 1 i ′ ) δ + (w 1 + w ′ 1 ) -(w 2 + w ′ 2 ) (3.18) = w 1 -w ′ 2 (3.19) so that P (i 2 = i ′ |(iii)) = 1 l {i ′ ∈R1} (k 2 i ′ ) δ w 1 + 1 l {i ′ ∈R2} (k 2 i ′ ) δ -(k 1 i ′ ) δ w 1 . (3.20)
(iv) If I ≥ w B + w 1 + w ′ 1 (this case cannot occur when δ = 1): then we do not define i 1 and i 2 .

Before going ahead with the definition of our coupling we note, as a direct consequence of our remarks in (i), (ii), (iii) and of of the fact that I has a uniform distribution: Proposition 3.1.1 For all i ∈ {1; . . . ; m}, P (i

1 = i|X 1 t , X 2 t ) = (k 1 i ) δ /l δ and P (i 2 = i|X 1 t , X 2 t ) = (k 2 i ) δ /l δ .
We then choose an integer j ∈ {1; . . . ; m} with uniform law and we distinguish once again between our four previous cases.

(i) If i 1 ∈ B: then i 2 = i 1 , we just write i 1 = i 2 = i. Then k 1 i = k 2 i . Thus, let a, b ∈ {1; 2} such that a = b and k a j ≤ k b j . Since both ψ [δ]
i and ψ

[δ] j are non-increasing,

p a := exp -ψ [δ] i (k a i -1) -ψ [δ] j (k a j ) + (3.21) ≥ exp -ψ [δ] i (k b i -1) -ψ [δ] j (k b j ) + =: p b . (3.22)
Let U be a uniform random variable on [0; 1).

-If U < p b : we set X a t+1 = (X a t ) ij and X b t+1 = (X b t ) ij . -If p b ≤ U < p a : we set X a t+1 = (X a t ) ij and X b t+1 = X b t . -If p a ≤ U : we set X a t+1 = X a t and X b t+1 = X b t .
In any of these cases, we certainly have ρ t+1 = ρ t .

(ii

) If i 1 ∈ R 2 : then i 2 = i 1 , we just write i 1 = i 2 = i. Let a, b ∈ {1; 2} such that a = b and ψ [δ] i (k a i -1) -ψ [δ] j (k a j ) ≤ ψ [δ] i (k b i -1) -ψ [δ] j (k b j ), so that p a := exp -ψ [δ] i (k a i -1) -ψ [δ] j (k a j ) + (3.23) ≥ exp -ψ [δ] i (k b i -1) -ψ [δ] j (k b j ) + =: p b . (3.24)
Let U be a uniform random variable on [0; 1).

-If U < p b : we set X a t+1 = (X a t ) ij and X b t+1 = (X b t ) ij . -If p b ≤ U < p a : we set X a t+1 = (X a t ) ij and X b t+1 = X b t . -If p a ≤ U : we set X a t+1 = X a t and X b t+1 = X b t .
In the last case we obviously have ρ t+1 = ρ t . In the first case the particles move together and ρ t+1 = ρ t . In the second case the number of red particles could increase only if

k a i ≤ k b i and k a j ≥ k b j , but, since ψ [δ]
i and ψ

[δ] j are non increasing, this would contradict p a > p b . As a consequence we have ρ t+1 ≤ ρ t in all the three cases.

(iii) If i 1 ∈ R 1 : there are three cases for i 2 . Either i 2 = i 1 = i and this case is the symmetric of (ii). Or i 2 is randomly chosen in R 2 , and we define

p 1 := exp -ψ [δ] i 1 (k 1 i 1 -1) -ψ [δ] j (k 1 j ) + , (3.25)
p 2 := exp -ψ [δ] i 2 (k 2 i 2 -1) -ψ [δ] j (k 2 j ) + . (3.26)
Or else i 2 is not defined, and we set p 2 := 0 still defining p 1 by (3.25). Let then V 1 , V 2 be independent uniform random variables on [0; 1).

-

If V 1 < p 1 : we set X 1 t+1 = (X 1 t ) i 1 j , else we set X 1 t+1 = X 1 t . -If V 2 < p 2 : we set X 2 t+1 = (X 2 t ) i 2 j , else we set X 2 t+1 = X 2 t .
In the first case we have ρ t+1 ≤ ρ t as previously. In the last two cases we also have ρ t+1 ≤ ρ t since only particles from R 1 in the first configuration and from R 2 in the second configuration can move.

(iv) If i 1 and i 2 are not defined: then we simply set (X 1 t+1 , X 2 t+1 ) = (X 1 t , X 2 t ) and we have ρ t+1 = ρ t .

In any of the previous cases, once i 1 , i 2 and j have been defined, the probability for X 1 t+1 (resp. X 2 t+1 ) to be (

X 1 t ) i 1 j (resp. (X 2 t ) i 2 j ) is exp{-[ψ [δ] i 1 (k 1 i 1 - 1) -ψ [δ] j (k 1 j )] + } (resp. exp{-[ψ [δ] i 2 (k 2 i 2 -1) -ψ [δ] j (k 2 j )] + }).
Thus, according to Proposition 3.1.1, the fact that j is uniformly chosen in {1; . . . ; m} and our study on the variation of ρ we conclude: Proposition 3.1.2 The process (X 1 , X 2 ) is a coupling for the dynamics defined by (1.35)-(1.36) and such that (ρ t ) t∈N is non-increasing.

Estimating the coupling time

Similarly to the proof of Theorem 1 we will use Proposition 2.2.1: since (ρ t ) t∈N is non-increasing it will be enough to give a lower bound for the probability of {ρ t+1 < ρ t }. Proposition 3.2.1 If at step t of the coupled dynamics (X 1 , X 2 ), we assume that "red particles have been chosen", i.e., i 1 ∈ R 1 and i 2 ∈ R 2 , then, there is a choice of j ∈ {i 1 ; i 2 } for which the number of red particles decreases with probability 1.

Proof: Assuming i 1 ∈ R 1 and i 2 ∈ R 2 yields k 1 i 1 > k 2 i 1 and k 1 i 2 < k 2 i 2 .
Using exactly the same argument as for Proposition 2.4.1 we prove that either exp{-[ψ

[δ] i 1 (k 1 i 1 -1) -ψ [δ] i 2 (k 1 i 2 )] + } = 1 or exp{-[ψ [δ] i 2 (k 2 i 2 -1) -ψ [δ] i 1 (k 2 i 1 )] + } = 1.
Eventually, if one red particle in some configuration moves to a site with a red particle in the other configuration, then both particles turn blue and the number of red particles decreases by one. Corollary 3.2.2 At step t, the probability for ρ t+1 to be ρ t -1 is at least δk δ-1 ρ t /ml δ .

Proof: The probability to choose i 1 ∈ R 1 and i 2 ∈ R 2 is

P (i 1 ∈ R 1 , i 2 ∈ R 2 ) = i ′ ∈R2 i∈R1 P (i 2 = i ′ |i 1 = i) × P (i 1 = i) (3.27) = i ′ ∈R2 i∈R1 (k 1 i ) δ -(k 2 i ) δ (k 1 i ) δ (k 2 i ′ ) δ -(k 1 i ′ ) δ λ(T ) × (k 1 i ) δ l δ (3.28) = 1 l δ i ′ ∈R2 ((k 2 i ′ ) δ -(k 1 i ′ ) δ ). ( 3
) ∈ R s + → f ( i z i ) -i f (z i
) is non-increasing in all its s variables (as a consequence of the slope inequalities), by concavity of z → z δ and using the fact that, for all i

′ ∈ R 2 , k 2 i ′ > k 1 i ′ we get P (i 1 ∈ R 1 , i 2 ∈ R 2 ) ≥ 1 l δ i ′ ∈R2 k 2 i ′ δ - i ′ ∈R2 k 1 i ′ δ (3.30) = 1 l δ ρ t + i ′ ∈R2 k 1 i ′ δ - i ′ ∈R2 k 1 i ′ δ . (3.31)
Using the same property of concave functions on R + (with s = 2) and the fact that i ′ ∈R2 k 1 i ′ ≤ k -ρ t , then using once again the concavity of z → z δ , we write

P (i 1 ∈ R 1 , i 2 ∈ R 2 ) ≥ 1 l δ k δ -(k -ρ t ) δ ≥ δk δ-1 ρ t l δ (3.32)
so that, by the previous proposition,

P (ρ t+1 = ρ t -1|X 1 t , X 2 t ) ≥ δk δ-1 ρ t ml δ . (3.33) As a consequence E(ρ t+1 |X 1 t , X 2 t ) ≤ 1 - δ km(k ∧ m)

A Simulation appendix

In this appendix we report on simulations for the temporal evolution of d(t) defined in (1.10) as well as the dependency in β of the mixing time t ǫ defined in (1.9) with ǫ = 0.1 in the context of the Fermi statistics ν associated with k = 50 particles among m = 20 energy levels with two different profiles.

• In the first case we chose v j = j/20 and n j = 2 j for 1 ≤ j ≤ m so that n = j 2 j = 2 21 -2. See Figure 1. Figure 1: Degeneracies of the energy levels for case 1, in semi-logarithmic representation for the second picture.

• In the second case we chose the same v j 's, the same volume n = 2 21 -2 and we chose the degeneracies according two a multinomial law with parameters n, 1/20, . . . , 1/20. See Figure 2. We were then interested in

d(t) = max η∈X k,m p t (η 0 , •) -ν T V = max η∈X k,m 1 2 η∈X k,m |p t (η 0 , η) -ν(η)| (A.1)
for t ≥ 0. The cardinality of X k,m being very large (equal to 69 19 ≃ 5 × 10 16 in case 2), one faces three difficulties with such a formula: will come back on this point later), we also considered a different η 0 , that one for which the high energy levels are empty and the low energy levels are saturated or as close to saturation as possible:

η 0 =
(2, 4, 8, 16, 20, 0, . . . , 0) in case 1, (50, 0, . . . , 0) in case 2. (A.7)

Then we simply took the maximum of the two quantities φ(µ N t ) -φ(ν N ) T V computed, in each case, with these two specific choices.

Summing up, we approximated d(t) by

dN 0 (t) = max η0∈X 0 k,m φ(µ N t ) -φ(ν N ) T V (A.8)
where X 0 k,m is the two configuration set described above. We chose N = 1024 and M = √ N = 32 and we plotted for different temperatures a graphical representation of the law of φ(ν N ) and the temporal evolution of dN 0 (t) for 0 ≤ t ≤ 2T . See Figures 3 to 7 For each temperature we get in each case an estimation tǫ of t ǫ . Estimating t ǫ in this way for 51 temperatures (β = 0 included) we show in Figure 13 the graphical representations of the simulated dependency of tǫ on β > 0 together with our upper bound from Theorem 1.

We can now conclude with a few comments. First, with such a procedure, t ǫ is essentially underestimated. Indeed, except for the replacement of p t (η 0 , •) and ν by µ N t and ν N our successive approximations underestimated d(t) then t ǫ . As a consequence we recover with case 2 that km ln(k/ǫ) is the best upper bound on t ǫ that is uniform over the disorder.

Second, we note that β can be seen as a tuning parameter for the disorder. We recover that our theoretical estimate on t ǫ can be improved by a factor of order m in weak disorder situations.

Last, it is worth to note that case 1 suggests that t ǫ is not generally decreasing with the temperature. In our preliminary simulations we only considered η 0 = (0, . . . , 0, 50) to estimate the mixing time. Then we feared that our non monotonic estimation could be an artefact of our simulation. This led us to introduce the alternative saturated low level starting configuration, confirming in that way the non-monotonicity.

1. 1

 1 Sampling the Fermi statisticsGiven two positive integers k and m, given a non-negative real number β, given m real numbers v 1 , . . . , v m and given m integers n 1 , . . . , n m such that

2 η∈X k,m |ν 1

 21 η∈X k,m p t (η, •) -ν TV (1.10)where • TV stands for the total variation distance defined for any probability measures ν 1 and ν 2 on X k,m byν 1 -ν 2 TV := max A⊂X k,m |ν 1 (A) -ν 2 (A)| = 1 (η) -ν 2 (η)|. (1.11)

1 l

 1 order to prove Theorem 1, we introduce a dynamics on the set of the possible distributions of the k particles in the m energy levels. Therefore we define the set of the configurations of k distinguishable particles in m energy levels Ω := {ω : {1; . . . ; k} → {1; . . . ; m}} (2.16) and for every such distribution ω, we define ξ(ω) = (k 1 , . . . , k m ) ∈ X k,m by ∀i ∈ {1; . . . ; m}, k i := k x=1 {ω(x)=i} .

Figure 2 :

 2 Figure 2: Degeneracies of the energy levels for case 2, with shifted and rescaled axis for the second picture.

Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 3456789101112 Figure 3: Probability distribution function of φ(ν N ) and temporal evolution of dN 0 (t) for case 1 and β = 0.

Figure 13 :

 13 Figure 13: Dependency of tǫ on β in case 1 and 2 respectively (semi-logarithmic representation). The horizontal line stands for our theoretical upper bound.

  for case 1 and Figures 8 to 12 for case 2.

Acknowledgment

J.R. is very grateful to everyone who made his stay in Rome possible and pleasant.

 [START_REF] Bakry | Diffusions hypercontractives, Séminaire de probabilités[END_REF]The whole research for this paper was done at the Dipartimento di matematica dell'università Roma Tre, during A.G.'s post-doc and J.R.'s research internship. This work was supported by the European Research Council through the "Advanced Grant" PTRELSS 228032.

i) One cannot compute ν(η) and p t (η 0 , η) for given η and η 0 .

ii) The sum in (A.1) contained too many terms to be computed.

iii) It is not possible to try all the possible η 0 before taking the maximum in (A.1).

We addressed the first difficulty by replacing p t (η 0 , •) by the empirical measure

computed with a large number N of simulated independent copies (X l t ) t∈N of our Markov chain starting from η 0 and, using Theorem 1, by replacing ν by

with T := ⌊km ln(k/ǫ)⌋ and where µ ′ t N and µ ′′ t N are independently computed like µ N t , starting from two different configurations. We addressed the second difficulty by looking at a coarse grained version Xk,m of X k,m through the free entropy φ. More precisely, with

we extended this partition of I N into a partition of R in intervals of length ∆ and grouped in a same class the configurations η ∈ X k,m that fall into a same interval. (For possible cases where φ -= φ + , we used a partition of R into intervals of length 10 -6 .) With φ the canonical projection from X k,m to Xk,m we could then compute φ(µ N t ) -φ(ν N ) T V rather than

and we also have, for large M and as soon as φ is an injection from X k,m to R,

As far as the third difficulty was concerned we underestimated

by making a guess on a particular configuration η 0 for which φ(µ N t )-φ(ν N ) T V could be of the order of dN (t). We simply took an η 0 for which we could have say, for all β, that it was "very far from typical configurations". We first took η 0 = (0, . . . , 0, k) since typical configurations are concentrated on the low energy levels when β is large and have their occupation numbers k j distributed like the degeneracies n j when β is small. Since this last point tends to show that in case 1 our specific η 0 is "not so far" from typical equilibrium configurations (we


Antoine Clément

Stéphane Laurens
email: stephane.laurens@insa-toulouse.fr

An alternative to Lyapunov Exponent as Damage Sensitive Feature

Keywords:

When using Lyapunov exponents as damage sensitive feature (DSF) in vibration based structural health monitoring, the computation time represents an important issue for real time applications. To overcome this problem, a new feature, the Jacobian Feature Vector (JFV), is proposed based on the same algorithm but allowing faster computation. Whereas for Lyapunov exponents calculation, the trajectories in the state space are followed for more than a thousand time steps, the new feature is formed by the Jacobian matrix of the dynamics after few time steps. The new feature and the Lyapunov exponents are tested on two case studies : a mass/spring/damper model and a laboratory three story structure. Both the DSF are showing a good performance in detecting damage for the numerical study, even if the Lyapunov exponents present a slightly better ability. Regarding the experimental data, Lyapunov exponents failed to detect damage while the JFV is clearly able to reveal structural changes. Therefore, the JFV offers a robust alternative to Lyapunov exponents allowing to save computation time.

Introduction

Structural Health Monitoring (SHM) consists in studying information which is generated by the structure to be investigated. The paradigm of non linear dynamical systems offers some interesting abilities for this purpose. From a general point of view, the scientific objective is focused on the detection of novelty in monitoring data. Among the various features available to describe dynamical systems, the Lyapunov Exponents (LE) are often proposed to detect changes in the system behavior.

A dynamical system can be represented by the evolution of its trajectory in its state space. Based on this representation, it is possible to calculate several invariants which characterize the dynamical system. Among them, Lyapunov exponents express the rate of divergence (or convergence) of trajectories initially close. If the system changes, due for example to an evolution of its parameters , the invariant property guaranties that the Lyapunov Exponents will deviate as well.

In the field of SHM, it is necessary to propose features that are sensitive to structural changes so that damage can be detected. The discrimination of damaged state is done by comparison with a reference baseline, generated under undamaged condition. The use of feature based on state space representation for time series analysis is limited in SHM. However, several studies explore the potential of this approach.

A first possibility is to excite the structure with a chaotic signal and to analyse the response in the state space. This method is used in several studies [START_REF] Nichols | Using state space predictive modeling with chaotic interrogation in detecting joint preload loss in a frame structure experiment[END_REF][START_REF] Todd | Vibration-based damage assessment utilizing state space geometry changes: local attractor variance ratio[END_REF] working on laboratory test structures, with different original features based on attractor geometry changes. Nichols et al [START_REF] Nichols | Structural health monitoring of offshore structures using ambient excitation[END_REF] show that it is possible to apply the same paradigm on a structure using ambient vibrations. A comparison of several metrics calculated on trajectory evolution for chaotic and broadband excitations has been carried out by Olson et al [START_REF] Olson | Sensitivity and computational comparison of statespace methods for structural health monitoring[END_REF].

The sensibility of Lyapunov exponents to damage has also been studied. They are used in the work reported by Ghafari et al [START_REF] Ghafari | Effect of localized faults on chaotic vibration of rolling element bearings[END_REF] in which the chaotic vibrations of rolling element bearings are investigated. Bearing failure is detected by the calculation of the largest Lyapunov Exponent.

Some studies have also proven the ability of Lyapunov Exponent to detect damage even in a non chaotic system. Livingston et al [START_REF] Livingston | Application of nonlinear dynamics analysis to damage detection and health monitoring of highway structures[END_REF] simulate cracks in a bridge FEM model excited by a traffic load model. The largest Lyapunov exponent is calculated along the bridge deck and is able to detect and locate the cracks. An application in a real structure have been carried out by Casciati et Casciati [START_REF] Casciati | Structural health monitoring by lyapunov exponents of nonlinear time series[END_REF] on a monumental arch (masonry specimen) before and after the repair of several major cracks. Two excitations are considered : hammer impact, and ambient vibrations. The method combines all the sensors to form the state space. It shows a good ability in detecting damage and even in locating cracks, as long as major non-stationarities are removed.

The last cited studies focus on the largest Lyapunov exponent or only the positives ones, but none of them uses the full spectrum. The main reasons are linked to the difficulty to estimate the negative Lyapunov exponents and to the high computation time. Though, when the aim is to find a Damage Sensitive Feature (DSF), even if the negative exponents estimations do not reflect exactly the real physical ones, they still may contain some information concerning the evolution of the system that can be useful in a comparative purpose.

This paper proposes a novel DSF based on Lyapunov exponent calculation algorithm : the Jacobian Feature Vector (JFV). It has to overcome the limitation of calculation time and has to minimize the loss of information about the dynamical system. After an introduction to the theoretical background leading to the definition of the JFV and LE, their practical calculation is presented. A comparison of the two DSF is then conducted on two case studies : 4 dof mass/spring/damper numerical model, and a laboratory structure. At last, a discussion on the performance of the DSF is proposed.

Theoretical background and DSF definition

The definition of the new DSF is closely related to the meaning of the LE.

Lyapunov exponents are linked to the predictability of the dynamical system studied. Viewed in the state space, let δ(0) << 1 be the distance between two nearby trajectories. After m time steps (m >> 1), the distance will be

δ(m) ≈ δ(0)e λm (1)
with λ the maximum Lyapunov exponent [START_REF] Kantz | Nonlinear Time Series Analysis[END_REF]. If the time is infinite, the ergodic theorem of Oseledec [START_REF] Oseledec | A multiplicative ergodic theorem: Lyapunov characteristic numbers for dynamical systems[END_REF] proves the existence and uniqueness of λ, whatever the initial point chosen. The Eq.1 illustrates the fact that if a small perturbation is imposed to the system, the rate of amplification of this perturbation is the Lyapunov exponent. If λ is positive, the trajectories will diverge making the system unpredictable after some time. This kind of behavior characterizes chaotic systems, which are sensitive to small perturbation. The calculation of the integer spectrum is done by following two initially close trajectories by the linearized dynamics [START_REF] Henry | The analysis of observed chaotic data in physical systems[END_REF].

Let X(k) ∈ R n be a point moving in a n-dimensional space at discrete time k = 1, 2, ..., N . The evolution of X(k) defines the trajectory and satisfies

X(k + 1) = F (X(k)) (2)
with F the dynamical system operator. If a nearby trajectory Y (k) ∈ R n is taken and noting δX i = X(k + i) -Y (k + i), the first-order Taylor expansion of F around X(k) leads to

δX 1 = X(k + 1) -Y (k + 1) = F (X(k)) -F (Y (k)) = J (X(k)) [X(k) -Y (k)] + o X(k) -Y (k) = J (X(k)) δX 0 + o δX 0 (3)
with J, the Jacobian matrix of F . When the trajectories are followed for L steps ahead, the product of L Jacobian matrices needs to be calculated and it is noted

J L (X k ) (Eq. 4). δX L = J(X k ) • J(X k+1 ) • ... • J(X k+L-1 ) • δX 0 = L-1 j=0 J(X k+j ) • δX 0 = J L (X k ) • δX 0 (4)
The Oseledec matrix, J L (X k ). J L (X k ) T , has eigenvalues noted Λ i , i = 1...n. Then the Lyapunov exponents are

λ i = lim L→∞ 1 2L ln(Λ i ) (5)
This equation (Eq. 5) implies that each trajectory needs to be followed for a very long time leading to :

• a very long computing time

• the hiding of every short irregularities of the trajectories.

These irregularities are the markers of damage induced, for example, by shocks in a loosened bolt, or by opening and closing of cracks. One is interested in detecting this local changes for SHM purpose. To highlight these irregularities, local Lyapunov exponents can be computed by following trajectories but only for a short time [START_REF] Abarbanel | Local lyapunov exponents computed from observed data[END_REF]. This approach has been more used to study the predictability of chaotic system than for SHM applications. But these irregularities should be visible directly on the Jacobian matrices since, in both cases (global or local), Lyapunov exponents represent the eigenvalues of the product of Jacobian matrices (Eq.5).

To preserve the potential sensitivity of local Lyapunov exponents and to provide a more simple approach, Jacobian matrix for short time evolution can be directly used to detect local distortions in trajectories. The components of the matrix are forming a new feature vector which will be referred to as the Jacobian Feature Vector (JFV).

Practical implementation

In experimental structural analysis, only a few degree of freedom of the system can be measured corresponding to each sensor location. Each measure is presented in the form of a time series. This section describes the algorithm applied on a single time series to obtain the Lyapunov spectrum and the JFV. The first step is to transform a scalar time series into a multi-dimensional time series to reconstruct the state space.

State Space Reconstruction

Based on the measure of an unique scalar time series (x) = {x(1), x(2), . . . , x(N )}, Takens proves that it is possible to qualitatively reconstruct the attractor in a ndimensional state space by using delayed coordinates method [START_REF] Takens | Detecting strange attractors in turbulence[END_REF]. This supposes that one single degree of freedom encloses enough information about the others.

Theoretically, the delay, τ , could be any entire number. In practice, it is common to choose τ long enough, to limit the redundancy between coordinates, and not too much, to keep some correlation. Once this choice is made, the state space vectors are formed, for k = 1, ..., n -τ (n -1), as follows:

X(k) = {x(k), x(k + τ ), ..., x(k + (n -1)τ )} (6)
Two methods for evaluating the optimal time delay τ are commonly used, the first minimum of the mutual information function and the first zero of the autocorrelation function [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF] . Then, one needs to find the proper embedding dimension, n, in which the attractor is completely unfolded. Brown et al suggest the method of false neighbors [START_REF] Brown | Computing the lyapunov spectrum of a dynamical system from an observed time series[END_REF] which counts the number of neighbors which do not remains neighbors as the dimension increases. The stabilization of this number indicates that the embedding dimension is large enough.

Once n and τ are selected, the time series is transformed into a collection of state space vectors according to Eq.6.

JFV and Lyapunov exponents Calculation

To evaluate the Lyapunov spectrum, the method described by Brown et al [START_REF] Brown | Computing the lyapunov spectrum of a dynamical system from an observed time series[END_REF] is used. The first step consists in calculating the Jacobian matrix along a trajectory. For this purpose, a fiducial point X(k) is chosen, and its r nearest neighbors , noted X r neigh (k), are selected. Difference vectors are formed according to Eq. 8.

X r neigh (k) = {X(i), i ∈ [1; N ] | X(k) -X(i) ≤ ∈ R} (7) δX r neigh (k) = X(k) -X(i) | X(i) ∈ X r neigh (k) (8)
To estimate the Taylor expansion of F , the best polynomial mapping, noted T (X), between each of the r perturbation vectors at time k and its evolutions at time k + T 2 (Eq.9) is calculated by a least-square method. Then, the Jacobian matrix associated to the linear part of T (X) is extracted.

δX i neigh (k + T 2 ) = T δX i neigh (k) (9)
with i = 0, . . . , r

The original approach uses a linear mapping, but the polynomial one is more suitable to estimate the negative exponents accurately [START_REF] Sano | Measurement of the lyapunov spectrum from a chaotic time series[END_REF]. The Jacobian matrix is estimated for LT 2 time steps, giving J(X k+i T 2 ), i = 1, ..., L.

Since the JVF has to reveal a local information, it is formed by the last Jacobian matrix for very few time steps. The L values from 1 to 5 time steps are tested (Eq.10).

JF V = J(X k+i T 2 )(:) i ≤ 5 (10)
On the other hand, to find the Lyapunov spectrum, the eigenvalues of the product of J(X k+i T 2 ) for i = 1, ..., L need to be calculated (Eq.4). L has to be at least 2000 to ensure the convergence of the limit (Eq.5). Since this product is very ill-conditioned, a QR decomposition method is recommended [START_REF] Eckmann | Liapunov exponents from time series[END_REF].

J(X k+(i+1)T 2 ).Q(i) = Q(i + 1).R(i + 1) (11)
with i = 0, . . . , L -1

Since Q(i) is an orthogonal matrix and R(i) an upper triangular one, the product of Eq. 4 is decomposed as

J L (X k ) = Q(L) L j=1 R(j) (12)
Finally, the Lyapunov exponents are expressed as

λ i = 1 dt L T 2 L j=L 0 ln (R ii (j)) i = 1, ..., n (13)
with dt the sampling period (in seconds) of the original scalar time series. The sum should start at L 0 = 0 to omit the first transient values.

For both Lyapunov exponents and JFV, the algorithm is repeated for several fiducial points across the attractor and the average of each DSF is taken.

SHM and Novelty Detection

For the long term monitoring of structures, damage can only be detected by comparison with a reference normal state acquired from various conditions [START_REF] Worden | The fundamental axioms of structural health monitoring[END_REF]. For each condition, time series are collected thanks to vibration measurements and compressed into feature vectors such JFV and LE. To compare a new feature vector to the baseline vectors collected at the reference state, the Mahalanobis distance is an appropriate approach [START_REF] Worden | Damage detection using outlier analysis[END_REF]. It gives the distance between a point and a set of points, taking into account the way the axis are populated. In statistics, it is commonly used to detect abnormal points in a data set. The advantage is that it provides a scalar feature from multivariate data. The Mahalanobis distance of a new feature vector DSF i (which can be either the Lyapunov exponents or the JVF) to a reference baseline characterized by a mean vector DSF and a covariance matrix C DSF is

D 2 i = DSF i -DSF T (C DSF ) -1 DSF i -DSF . ( 14
)
To quantify the efficiency of the features proposed for damage detection, the number of misclassified points is calculated. The database, composed of time series collected at the undamaged state, and time series collected for various damage levels, is split in three parts :

• the reference baseline, composed of half of the undamaged database. It is used as reference to form the covariance matrix (C DSF ) and mean vector (DSF ) for Mahalanobis distance calculation. • the test base, composed of the last quarter of the undamaged database and all the damaged time series. It is used to perform a blind-test aimed at checking the efficiency of each DSF. In this study, since each vector is known, the efficiency of the DSF will be quantified by the number of misclassified vectors, including :

undamaged state classified as damaged (False Alarm) damaged state classified as undamaged (Non-Detection)

Case Studies

Case Study 1

The first case study (CS1) is a 4 dof mass/spring/damper numerical model (Fig. 1). The physic constants are selected to be realistic with building structures : k = 2.4e8 N.m -1 , c = 1.02e6 N.m -1 .s and m = 70000 kg. The system can be expressed by the equation

M ẍ + C ẋ + Kx = f(t) (15)
where M , C and K represent mass, damping and stiffness 4 × 4 matrices, respectively. f(t) is the excitation vector which applies a white noise at the free end of the model. Equation 15 is solved by a fourth order Runge-Kutta method sampled at 60Hz. The displacement x 1 , resampled at 30Hz, is used for damage detection. An example of time series calculated is visible in Fig. 2.

Damage is simulated by a decrease in the first spring stiffness, but only when it is in extension (Eq.16) to represent the effect of the closing and opening of a crack.

k 1 = k if x 1 ≥ 0 (1 -α)k if x 1 < 0 ( 16
)
As presented in Tab.1, the database is composed of 120 undamaged time series and 90 damaged time series with an increasing level of damage.

Case Study 2

For the second case study (CS2), vibration signals of a three story bookshelf structure are considered (Fig. 3). These experiments are a part of several case studies carried out by the Los Alamos National Laboratory (LANL) and provided freely for the scientific community. All the informations about the structure and the test set up are available on LANL website [START_REF]The bookshelf structure[END_REF].

It is composed of Unitrust columns and aluminum floor plates. Support brackets for the columns are bolted to the plates and hold the Unitrust columns. The global geometry is detailed in Fig. 4.

A shaker is attached at corner D (Fig. 4) in order to excite both torsional and translational vibration modes. Twenty-four accelerometers are positioned on the structure, all in the same direction. The excitation is white noise with several bandwidths from 800 to 3200 Hz, and different input levels from 3,5 to 7 Volts. For each setup, the response is sampled at 1600 Hz and consists of 8192 points.

For the purpose of testing the sensitivity to damage, only one sensor is considered. It is located at the corner C, on the first floor plate. This is close to the location of the induced damage.

Damage is simulated by changes in the corner C first floor joint. Tab.2 summarizes the different damage cases. The database consists of 150 time series related to the undamaged state, and 60 time series in the damaged state.

Parameters investigated

In the algorithm described previously, several parameters need to be specified. The first ones are the reconstruction parameters. They are set based on the analysis of baseline series chosen randomly. The first minimum of the mutual information function and the first zero of the autocorrelation function indicate that a delay of 2 for case study 1, and 6 for case study 2, are good choices. However, if T 2 is equal to τ , the two vectors mapped by the polynomial fitting are almost identical leading to identity matrix as Jacobian. Since the maximum value investigated is T 2 = 5, the choice of the delay is imposed by T 2 for CS1. Finally, τ is set to 6 for both case studies. The false neighbors method indicates that a dimension of 7 for CS1 and CS2 is a good compromise between the stabilization of the number of false neighbors and a small dimension which keeps the calculation time low (Fig. 5).

It is important to note that there is no need to have the same embedding parameters for the two case studies. For Lyapunov spectrum determination, a preliminary study has shown that the number of transient mapping steps does not exceed 500 and that they are stable after 1500 other ones. These results lead to the choice of L = 2000 and L 0 = 500. Each neighborhood is composed of N para to 2N para points, with N para the number of parameters of the polynomial mapping function to be fitted (depending on the order). At last, 40 fiducial points are followed in each attractor for Lyapunov exponents calculation. For the JFV, since the averaging over the whole attractor is not done according to a large L, the number of fiducial points will be 100.

The order of the polynomial mapping and the mapping steps (T 2 ) are supposed to be influencing parameters. Therefore, several values are tested. However, the order have to be smaller than 4 to limit the number of points in each neighborhood (329 points for order 4). A similar problem appears with T 2 in Lyapunov exponents determination, it can not be higher than 3 since the time series are composed of 8192 points and the trajectories are followed for 2000 mapping steps. In the case of the JFV, the number of mapping steps before extracting the Jacobian matrix is investigated. Tab. 3 summarizes the values of the parameters involved in the study, leading to 16 parameter setups.

Results and discussion

Lyapunov Exponents

The efficiency of Lyapunov Exponents as DSF is very different between the two case studies. For CS2, the percentage of misclassified points varies from 41% to 62% (Tab.4) including no false alarms, but a large amont of non-detection (Fig. 6b). On the other hand, for CS1, LE are capable of detecting damage from 30%-40% of stiffness reduction, depending on the set of parameters. The poor results obtained in CS2 is caused by the dispersion of the LE values observed on experimental data. This dispersion makes difficult the separation between damaged and undamaged states (Fig. 7b).

To determine accurate Lyapunov exponents of chaotic systems, the use of a high order polynomial mapping is preferred. But Tab.4 shows that increasing the order of the mapping leads to more misclassified points. When the polynomial order increases, it is observed that the variations of the different exponents tend to be uniform, leading to a decrease in the sensitivity since each exponent carries the same information (Fig. 7). Higher order also produces several positive exponents. Since Lyapunov exponents of noise are theoretically infinite, it probably means that the random part of the time series is modeled by the mapping. The results are not improved by doubling the mapping step.

Jacobian Feature Vector

The number of misclassified time series for the different sets of parameters is presented in Tab.5. Even if some observations will be formulated regarding the influence of the parameter setup, the classification ability is comparable from one setup to another.

For CS1, the worst case counts 51% of misclassified points and 21% for the best one, whereas for CS2 the worst case is 16% which is far better. The change in spring stiffness is not detected before 40% of reduction, whereas all damage cases are correctly classified for CS2 (Fig. 8). This difference is explained by the progressivity of damage introduction in CS1 which make the early detection more challenging.

Looking at the results in details reveals some trends in the parameters influence. Tab.5 shows that increasing the order does not reduce the number of misclassified points. On contrary, it increases for CS2, even if it remains small. The explanation is probably the same as for the Lyapunov exponents. For damage detection, the first order appears more relevant.

The number of misclassified points as T 2 increases does not present a clear trend. When L = 1 (setups #4 -5 -6), large T 2 improves the sensitivity for CS1, but not for CS2. Whereas with L = 2 (setups #7 to 11), results remain stable, except for T 2 = 3 (setup #9) which will be explained later. This proves that mapping step has no direct effect on the sensitivity to damage, and that it has to be considered in interaction with L.

Change in the value of L, implies irregular variations in the efficiency of the JFV as presented in Tab.5. For both CS1 and CS2, there are no trends as L increases. However In order to understand why, for these particular cases, the sensitivity of the JFV is reduced, one has to consider how the different vectors of the state space are reconstructed. The optimal reconstruction delay selected is six time step (τ = 6). The Fig. 9 shows the ten first vectors reconstructed from the time series (x) = {x 1 , x 2 , ..., x N }. In the case T 2 = 2, L = 4, the Jacobian matrix used in the JFV is calculated between the vectors X 7 and X 9 which differ only by one coordinate from initial vector X 1 and vector X 3 . Therefore this case is equivalent to the case L = 1 which presents a reduced sensibility. By the same token, the case T 2 = 2, L = 3 uses a mapping between X 5 and X 7 . It is close to a mapping between X 1 and X 5 formulated with L = 1 and T 2 = 4. The same explanation holds for the case T 2 = 3 and L = 2. The Jacobian matrix of the first mapping step (L = 1) is less sensitive to damage than the other ones because the initial neighbors are much more close than after few time steps. If the points are spread in a larger area in the attractor, the Jacobian matrix is more able to capture changes in the dynamics.

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 10 x 1 x 2 x 3 x 4 x

Discussion

The sensitivity of the JFV depends on the set of parameters chosen. As explained previously, for an optimal sensibility, it is important to avoid the inclusion of the initial neighborhood (and its reconstruction multiples) in the Jacobian calculation for the JFV. This condition is verified if

T 2 • L = τ and (17)
T 2 • L = τ -L with L > 1 (18)
Keeping T 2 as small as possible will contribute to reduce the processing time. In this study, this is achieved with (Fig. 8a et Fig. 8c) :

• linear mapping

• trajectories followed for two mapping steps (L = 2)

• mapping steps formed by one or two time steps (T 2 = 1, 2).

Once these conditions are verified, the JFV is able to detect damage quite effectively. For the CS1, in which data are obtained with a numerical model, the efficiency of the LE is slightly better than the JFV since 30% of stiffness reduction is detected by LE, against 40% by the JFV. However, on real data, LE hardly detect the damage introduced on the bookshelf structure, whereas all cases are correctly identified by JFV. It is likely that time series in CS2 contain some noise due to experimental conditions which is not present in the simulated data. It seems that JFV is more robust to noise than LE. Moreover, to calculate the Lyapunov spectrum of a time series with an Intel R Xeon 3GHz processor, the CPU time varies from 170s to 3700s depending the polynomial order, whereas it is only between 0.2s to 6s for the new JFV.

The study shows that the proposed feature is an interesting alternative to Lyapunov exponents for SHM purpose, regarding both sensitivity to damage and computation time.

Conclusion

Using the full Lyapunov spectrum as damage sensitive feature is very time consuming. To overcome this issue, the Jacobian Feature Vector (JFV), a new damage sensitive feature based on the Jacobian matrix of reconstructed dynamics in the state space, is proposed. Instead of following the trajectories for several thousand time steps, as needed for Lyapunov Exponents (LE) calculation, the JFV involves the Jacobian matrix at the early time steps. Both JFV and LE are compared to detect damage on a simulated mass/spring model and a laboratory bookshelf structure. Unlike LE, the JVF is able to detect the damage introduced in the two case studies with an appropriate set of parameters. This indicates that the proposed approach is more robust than LE. Future investigations will be focused on :

• the verification of the sensibility to damage with variable signal-to-noise ratios

• other case studies

• the comparison of the JFV approach to common SHM paradigms (modal analysis, autoregressive models etc.).

Figure 1 .

 1 Figure 1. Mass/spring/damper model

Figure 2 .

 2 Figure 2. Time series obtained with the model.

Figure 3 .

 3 Figure3. Photo of the test structure (LANL)[START_REF]The bookshelf structure[END_REF]

Figure 4 .

 4 Figure 4. Front (left figure) and top dimensions (right figure) of the test structure (LANL)[START_REF]The bookshelf structure[END_REF]

Figure 5 .

 5 Figure 5. Baseline series #5 : Number of false neighbors functions vs embedding dimension both case studies

Figure 6 .

 6 Figure 6. Lyapunov Exponents Mahalanobis distance of the test base for T 2 = 1 and linear mapping. The dashed line represents the classification threshold. a) CS1 ; b) CS2 ;

Figure 7 . 2 •

 72 Figure 7. Lyapunov spectrum for the undamaged time series (#1 to #150 for CS1 and #1 to #120 for CS2) and the damaged time series. Left figure display LE for CS1 and right for CS2. a1) and b1) order= 1 ; a2) and b2) order= 2 ; a3) and b3) order= 3 ;

Figure 8 .

 8 Figure 8. JFV Mahalanobis distance of the test base for the worst and best classification results. The dashed line represents the classification threshold. a) setup #1 CS1 Best ; b) setup #1 CS2 Best ; c) setup #4 CS1 Worst ; d) setup #4 CS2 Worst.

Figure 9 .

 9 Figure 9. Ten first reconstructed vectors.

Table 1 .

 1 Composition of the database

	Number of time series	α
	120	0
	10	0.1
	10	0.2
	10	0.3
	10	0.4
	10	0.5
	10	0.6
	10	0.7
	10	0.8
	10	0.9

Table 2 .

 2 Damage cases : modification of the bolts between the bracket and the plate (LANL)[START_REF]The bookshelf structure[END_REF]

	Damage Case Description	Bolt Torque [N.m]	Number of time series collected
	0	Undamaged	6.77	150
	1	Decrease torque	1.128	15
	2	Decrease tightening torque	0.564	15
	3	hand tighten	-	10
	4	Bolts removed	-	10
	5	Bracket removed	-	10

Table 3 .

 3 Parameters values for DSF calculation

	Feature	T 2	L	Order
	Lyapunov exponents	1 2	2000 2000	1,2,3 1,2
		2	2	1,2,3
	JFV	1,2,3 1,2,3,4,5 2 1	1 1
		2	1,2,3,4,5 1

Table 4 .

 4 Parameters values and percentage of misclassified points in the test-base for Lyapunov exponents

	T 2 Order % of Misclassified points
			CS1	CS2
	1	1	15	41
	1	2	23	53
	1	3	24	62
	2	1	25	61
	2	2	31	60

Table 5 .

 5 Parameters values and percentage of misclassified points in the test-base for JFV

	Set up number T 2 L Order % Misclassified points
					CS1	CS2
	1)			1	27	0
	2)	2	2	2	34	1
	3)			3	34	5
	4)	1			51	14
	5)	2	1	1	43	13
	6)	3			36	16
	7)	1			30	0
	8)	2			25	0
	9)	3	2	1	42	9
	10)	4			28	1
	11)	5			46	1
	12)		1		46	13
	13)		2		32	0
	14)	2	3	1	49	4
	15)		4		30	11
	16)		5		21	0

  5 x 6 x 7 x 8 x 9 x 10 x 7 x 8 x 9 x 10 x 11 x 12 x 13 x 14 x 15 x 16 x 13 x 14 x 15 x 16 x 17 x 18 x 19 x 20 x 21 x 22 x 19 x 20 x 21 x 22 x 23 x 24 x 25 x 26 x 27 x 29 x 25 x 26 x 27 x 28 x 29 x 30 x 31 x 32 x 33 x 34 x 31 x 32 x 33 x 34 x 35 x 36 x 37 x 38 x 39 x 40 x 37 x 38 x 39 x 40 x 41 x 42 x 43 x 44 x 45 x 46

Acknowledgments

The authors gratefully acknowledge the Los Alamos National Laboratory for the data provided to the scientific community and the Ministère de l'Enseignement Supérieur et de la Recherche for its funding.